Setting up the R workspace
library(tidyverse)
── Attaching core tidyverse packages ─────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ purrr 1.0.2 ✔ tibble 3.2.1
✔ readr 2.1.5 ✔ tidyr 1.3.1
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(rjags)
Loading required package: coda
Linked to JAGS 4.3.2
Loaded modules: basemod,bugs
library(ggplot2)
library(hydroGOF)
set.seed(1)
Read in the data here, which is an output of XX code file
Then separate the data into air and water temperature files
# Air and water temperature data: CR loccation, 4.5 scenario
CR_Temps <- read.csv("CR_45SNAPForecast.csv", header=TRUE)
# Air temperature variables
CR_AirTemp45 <- CR_Temps %>%
select(date, Airport, Air_MonthAvg, std, precision, Region)
# Water temperature variables
CR_WaterTemp45 <- CR_Temps %>%
select(date, BVS, CAB, TIN, SQR, CAN, EYS, TIS, WDD, RHM)
# Set up a string of all the dates
time <- CR_WaterTemp45$date
length(time)
[1] 1164
Setting up the random walk model
dlm_pooled <- "
model{
#### Priors
for(p in 1:np){
x[1,p] ~ dnorm(x_ic, tau_ic) # Initial condition of water temperature
}
tau_obs ~ dgamma(a_obs, r_obs) # Prior on observation error
tau_add ~ dgamma(a_add, r_add) # Prior on process error
#### Fixed Effects
beta ~ dmnorm(mu_beta, tau_beta) # Prior on beta coefficients
for(p in 1:np){
int[p] ~ dnorm(mu_int, tau_int) # Prior on pond-specific intercepts
}
#### Data Model
for(t in 1:n){ # loop over all time steps
for(p in 1:np){
OBS[t,p] ~ dnorm(x[t,p], tau_obs) # Observed water temperature is drawn from latent air temperature with observation uncertainty
}
Xf[t] ~ dnorm(muAirTemp[t], tauAirTemp[t]) # Latent air temperature is drawn from mean and precision of forecasated air temperature
}
#### Process Model
for(t in 2:n){ # loop over all time steps except teh first (we defined ic above)
for(p in 1:np){
mu[t,p] <- x[t-1,p] + int[p] + beta[1] * x[t-1,p] + beta[2] * Xf[t] # Mean water temperature is a function of the previous time step and current air temperature
x[t,p] ~ dnorm(mu[t,p], tau_add) # Latent water temperature is drawn from mean water temperature with process uncertainty
}
}
}
"
Define the data and priors for the model
# Empty list
data <- list()
# Water temperature observations
data$OBS <- dplyr::select(CR_WaterTemp45, -date)
# Number of time steps
data$n <- nrow(data$OBS)
# Number of ponds
data$np <- ncol(data$OBS)
# Initial water temperature mean
data$x_ic <- 0.1
# Initial water temperature precision
data$tau_ic = 0.1
# Prior parameters for observation and process uncertainty
data$a_obs = 1
data$r_obs = 1
data$a_add = 1
data$r_add = 1
# Prior parameters for beta coefficients
data$mu_beta <- c(0, 0)
data$tau_beta <- diag(x = c(0.001, 0.001), nrow = 2, ncol = 2)
# Prior parameters for intercept
data$mu_int <- 0
data$tau_int <- 0.001
# Mean air temperature estimate
data$muAirTemp <- CR_AirTemp45$Air_MonthAvg
# Air temperature precision
data$tauAirTemp <- CR_AirTemp45$precision
Create JAGS model with 3 chains
jm <- jags.model(file = textConnection(dlm_pooled), data = data, n.chains = 3)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 448
Unobserved stochastic nodes: 21680
Total graph size: 46573
Initializing model
Posterior samples of parameters
CRD45_out_params <- coda.samples(model = jm,
variable.names = c('beta', 'int',
'tau_add', 'tau_obs'),
n.iter = 50000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Posterior samples of response variables
CRD45_out_response <- coda.samples(model = jm,
variable.names = c('x', 'OBS'),
n.iter = 100000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
A couple quick checks here
# plot(CRD45_out_params)
gelman.diag(CRD45_out_params, confidence = 0.99)
Potential scale reduction factors:
Point est. Upper C.I.
beta[1] 1.00 1.01
beta[2] 1.01 1.04
int[1] 1.01 1.03
int[2] 1.01 1.04
int[3] 1.01 1.02
int[4] 1.00 1.02
int[5] 1.01 1.04
int[6] 1.01 1.04
int[7] 1.01 1.04
int[8] 1.01 1.04
int[9] 1.01 1.04
tau_add 1.00 1.00
tau_obs 1.00 1.00
Multivariate psrf
1.01
Visualize the output by just looking at the 95% credible interval of the time-series of X’s and compare that to the observed Y’s
Transform the samples back from the log domain to the linear domain
time ## adjust to zoom in and out
[1] "2012-01-01" "2012-02-01" "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01"
[8] "2012-08-01" "2012-09-01" "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01"
[15] "2013-03-01" "2013-04-01" "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01"
[22] "2013-10-01" "2013-11-01" "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01"
[29] "2014-05-01" "2014-06-01" "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01"
[36] "2014-12-01" "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01"
[43] "2015-07-01" "2015-08-01" "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01"
[50] "2016-02-01" "2016-03-01" "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01"
[57] "2016-09-01" "2016-10-01" "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01"
[64] "2017-04-01" "2017-05-01" "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01"
[71] "2017-11-01" "2017-12-01" "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01"
[78] "2018-06-01" "2018-07-01" "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01"
[85] "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01"
[92] "2019-08-01" "2019-09-01" "2019-10-01" "2019-11-01" "2019-12-01" "2012-01-01" "2012-02-01"
[99] "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01" "2012-08-01" "2012-09-01"
[106] "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01" "2013-03-01" "2013-04-01"
[113] "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01" "2013-10-01" "2013-11-01"
[120] "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01" "2014-05-01" "2014-06-01"
[127] "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01" "2014-12-01" "2015-01-01"
[134] "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01" "2015-07-01" "2015-08-01"
[141] "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01" "2016-02-01" "2016-03-01"
[148] "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01"
[155] "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01" "2017-04-01" "2017-05-01"
[162] "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01" "2017-11-01" "2017-12-01"
[169] "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01" "2018-06-01" "2018-07-01"
[176] "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01" "2019-01-01" "2019-02-01"
[183] "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01" "2019-08-01" "2019-09-01"
[190] "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01" "2020-03-01" "2020-04-01"
[197] "2020-05-01" "2020-06-01" "2020-07-01" "2020-08-01" "2020-09-01" "2020-10-01" "2020-11-01"
[204] "2020-12-01" "2021-01-01" "2021-02-01" "2021-03-01" "2021-04-01" "2021-05-01" "2021-06-01"
[211] "2021-07-01" "2021-08-01" "2021-09-01" "2021-10-01" "2021-11-01" "2021-12-01" "2022-01-01"
[218] "2022-02-01" "2022-03-01" "2022-04-01" "2022-05-01" "2022-06-01" "2022-07-01" "2022-08-01"
[225] "2022-09-01" "2022-10-01" "2022-11-01" "2022-12-01" "2023-01-01" "2023-02-01" "2023-03-01"
[232] "2023-04-01" "2023-05-01" "2023-06-01" "2023-07-01" "2023-08-01" "2023-09-01" "2023-10-01"
[239] "2023-11-01" "2023-12-01" "2024-01-01" "2024-02-01" "2024-03-01" "2024-04-01" "2024-05-01"
[246] "2024-06-01" "2024-07-01" "2024-08-01" "2024-09-01" "2024-10-01" "2024-11-01" "2024-12-01"
[253] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01" "2025-06-01" "2025-07-01"
[260] "2025-08-01" "2025-09-01" "2025-10-01" "2025-11-01" "2025-12-01" "2026-01-01" "2026-02-01"
[267] "2026-03-01" "2026-04-01" "2026-05-01" "2026-06-01" "2026-07-01" "2026-08-01" "2026-09-01"
[274] "2026-10-01" "2026-11-01" "2026-12-01" "2027-01-01" "2027-02-01" "2027-03-01" "2027-04-01"
[281] "2027-05-01" "2027-06-01" "2027-07-01" "2027-08-01" "2027-09-01" "2027-10-01" "2027-11-01"
[288] "2027-12-01" "2028-01-01" "2028-02-01" "2028-03-01" "2028-04-01" "2028-05-01" "2028-06-01"
[295] "2028-07-01" "2028-08-01" "2028-09-01" "2028-10-01" "2028-11-01" "2028-12-01" "2029-01-01"
[302] "2029-02-01" "2029-03-01" "2029-04-01" "2029-05-01" "2029-06-01" "2029-07-01" "2029-08-01"
[309] "2029-09-01" "2029-10-01" "2029-11-01" "2029-12-01" "2030-01-01" "2030-02-01" "2030-03-01"
[316] "2030-04-01" "2030-05-01" "2030-06-01" "2030-07-01" "2030-08-01" "2030-09-01" "2030-10-01"
[323] "2030-11-01" "2030-12-01" "2031-01-01" "2031-02-01" "2031-03-01" "2031-04-01" "2031-05-01"
[330] "2031-06-01" "2031-07-01" "2031-08-01" "2031-09-01" "2031-10-01" "2031-11-01" "2031-12-01"
[337] "2032-01-01" "2032-02-01" "2032-03-01" "2032-04-01" "2032-05-01" "2032-06-01" "2032-07-01"
[344] "2032-08-01" "2032-09-01" "2032-10-01" "2032-11-01" "2032-12-01" "2033-01-01" "2033-02-01"
[351] "2033-03-01" "2033-04-01" "2033-05-01" "2033-06-01" "2033-07-01" "2033-08-01" "2033-09-01"
[358] "2033-10-01" "2033-11-01" "2033-12-01" "2034-01-01" "2034-02-01" "2034-03-01" "2034-04-01"
[365] "2034-05-01" "2034-06-01" "2034-07-01" "2034-08-01" "2034-09-01" "2034-10-01" "2034-11-01"
[372] "2034-12-01" "2035-01-01" "2035-02-01" "2035-03-01" "2035-04-01" "2035-05-01" "2035-06-01"
[379] "2035-07-01" "2035-08-01" "2035-09-01" "2035-10-01" "2035-11-01" "2035-12-01" "2036-01-01"
[386] "2036-02-01" "2036-03-01" "2036-04-01" "2036-05-01" "2036-06-01" "2036-07-01" "2036-08-01"
[393] "2036-09-01" "2036-10-01" "2036-11-01" "2036-12-01" "2037-01-01" "2037-02-01" "2037-03-01"
[400] "2037-04-01" "2037-05-01" "2037-06-01" "2037-07-01" "2037-08-01" "2037-09-01" "2037-10-01"
[407] "2037-11-01" "2037-12-01" "2038-01-01" "2038-02-01" "2038-03-01" "2038-04-01" "2038-05-01"
[414] "2038-06-01" "2038-07-01" "2038-08-01" "2038-09-01" "2038-10-01" "2038-11-01" "2038-12-01"
[421] "2039-01-01" "2039-02-01" "2039-03-01" "2039-04-01" "2039-05-01" "2039-06-01" "2039-07-01"
[428] "2039-08-01" "2039-09-01" "2039-10-01" "2039-11-01" "2039-12-01" "2040-01-01" "2040-02-01"
[435] "2040-03-01" "2040-04-01" "2040-05-01" "2040-06-01" "2040-07-01" "2040-08-01" "2040-09-01"
[442] "2040-10-01" "2040-11-01" "2040-12-01" "2041-01-01" "2041-02-01" "2041-03-01" "2041-04-01"
[449] "2041-05-01" "2041-06-01" "2041-07-01" "2041-08-01" "2041-09-01" "2041-10-01" "2041-11-01"
[456] "2041-12-01" "2042-01-01" "2042-02-01" "2042-03-01" "2042-04-01" "2042-05-01" "2042-06-01"
[463] "2042-07-01" "2042-08-01" "2042-09-01" "2042-10-01" "2042-11-01" "2042-12-01" "2043-01-01"
[470] "2043-02-01" "2043-03-01" "2043-04-01" "2043-05-01" "2043-06-01" "2043-07-01" "2043-08-01"
[477] "2043-09-01" "2043-10-01" "2043-11-01" "2043-12-01" "2044-01-01" "2044-02-01" "2044-03-01"
[484] "2044-04-01" "2044-05-01" "2044-06-01" "2044-07-01" "2044-08-01" "2044-09-01" "2044-10-01"
[491] "2044-11-01" "2044-12-01" "2045-01-01" "2045-02-01" "2045-03-01" "2045-04-01" "2045-05-01"
[498] "2045-06-01" "2045-07-01" "2045-08-01" "2045-09-01" "2045-10-01" "2045-11-01" "2045-12-01"
[505] "2046-01-01" "2046-02-01" "2046-03-01" "2046-04-01" "2046-05-01" "2046-06-01" "2046-07-01"
[512] "2046-08-01" "2046-09-01" "2046-10-01" "2046-11-01" "2046-12-01" "2047-01-01" "2047-02-01"
[519] "2047-03-01" "2047-04-01" "2047-05-01" "2047-06-01" "2047-07-01" "2047-08-01" "2047-09-01"
[526] "2047-10-01" "2047-11-01" "2047-12-01" "2048-01-01" "2048-02-01" "2048-03-01" "2048-04-01"
[533] "2048-05-01" "2048-06-01" "2048-07-01" "2048-08-01" "2048-09-01" "2048-10-01" "2048-11-01"
[540] "2048-12-01" "2049-01-01" "2049-02-01" "2049-03-01" "2049-04-01" "2049-05-01" "2049-06-01"
[547] "2049-07-01" "2049-08-01" "2049-09-01" "2049-10-01" "2049-11-01" "2049-12-01" "2050-01-01"
[554] "2050-02-01" "2050-03-01" "2050-04-01" "2050-05-01" "2050-06-01" "2050-07-01" "2050-08-01"
[561] "2050-09-01" "2050-10-01" "2050-11-01" "2050-12-01" "2051-01-01" "2051-02-01" "2051-03-01"
[568] "2051-04-01" "2051-05-01" "2051-06-01" "2051-07-01" "2051-08-01" "2051-09-01" "2051-10-01"
[575] "2051-11-01" "2051-12-01" "2052-01-01" "2052-02-01" "2052-03-01" "2052-04-01" "2052-05-01"
[582] "2052-06-01" "2052-07-01" "2052-08-01" "2052-09-01" "2052-10-01" "2052-11-01" "2052-12-01"
[589] "2053-01-01" "2053-02-01" "2053-03-01" "2053-04-01" "2053-05-01" "2053-06-01" "2053-07-01"
[596] "2053-08-01" "2053-09-01" "2053-10-01" "2053-11-01" "2053-12-01" "2054-01-01" "2054-02-01"
[603] "2054-03-01" "2054-04-01" "2054-05-01" "2054-06-01" "2054-07-01" "2054-08-01" "2054-09-01"
[610] "2054-10-01" "2054-11-01" "2054-12-01" "2055-01-01" "2055-02-01" "2055-03-01" "2055-04-01"
[617] "2055-05-01" "2055-06-01" "2055-07-01" "2055-08-01" "2055-09-01" "2055-10-01" "2055-11-01"
[624] "2055-12-01" "2056-01-01" "2056-02-01" "2056-03-01" "2056-04-01" "2056-05-01" "2056-06-01"
[631] "2056-07-01" "2056-08-01" "2056-09-01" "2056-10-01" "2056-11-01" "2056-12-01" "2057-01-01"
[638] "2057-02-01" "2057-03-01" "2057-04-01" "2057-05-01" "2057-06-01" "2057-07-01" "2057-08-01"
[645] "2057-09-01" "2057-10-01" "2057-11-01" "2057-12-01" "2058-01-01" "2058-02-01" "2058-03-01"
[652] "2058-04-01" "2058-05-01" "2058-06-01" "2058-07-01" "2058-08-01" "2058-09-01" "2058-10-01"
[659] "2058-11-01" "2058-12-01" "2059-01-01" "2059-02-01" "2059-03-01" "2059-04-01" "2059-05-01"
[666] "2059-06-01" "2059-07-01" "2059-08-01" "2059-09-01" "2059-10-01" "2059-11-01" "2059-12-01"
[673] "2060-01-01" "2060-02-01" "2060-03-01" "2060-04-01" "2060-05-01" "2060-06-01" "2060-07-01"
[680] "2060-08-01" "2060-09-01" "2060-10-01" "2060-11-01" "2060-12-01" "2061-01-01" "2061-02-01"
[687] "2061-03-01" "2061-04-01" "2061-05-01" "2061-06-01" "2061-07-01" "2061-08-01" "2061-09-01"
[694] "2061-10-01" "2061-11-01" "2061-12-01" "2062-01-01" "2062-02-01" "2062-03-01" "2062-04-01"
[701] "2062-05-01" "2062-06-01" "2062-07-01" "2062-08-01" "2062-09-01" "2062-10-01" "2062-11-01"
[708] "2062-12-01" "2063-01-01" "2063-02-01" "2063-03-01" "2063-04-01" "2063-05-01" "2063-06-01"
[715] "2063-07-01" "2063-08-01" "2063-09-01" "2063-10-01" "2063-11-01" "2063-12-01" "2064-01-01"
[722] "2064-02-01" "2064-03-01" "2064-04-01" "2064-05-01" "2064-06-01" "2064-07-01" "2064-08-01"
[729] "2064-09-01" "2064-10-01" "2064-11-01" "2064-12-01" "2065-01-01" "2065-02-01" "2065-03-01"
[736] "2065-04-01" "2065-05-01" "2065-06-01" "2065-07-01" "2065-08-01" "2065-09-01" "2065-10-01"
[743] "2065-11-01" "2065-12-01" "2066-01-01" "2066-02-01" "2066-03-01" "2066-04-01" "2066-05-01"
[750] "2066-06-01" "2066-07-01" "2066-08-01" "2066-09-01" "2066-10-01" "2066-11-01" "2066-12-01"
[757] "2067-01-01" "2067-02-01" "2067-03-01" "2067-04-01" "2067-05-01" "2067-06-01" "2067-07-01"
[764] "2067-08-01" "2067-09-01" "2067-10-01" "2067-11-01" "2067-12-01" "2068-01-01" "2068-02-01"
[771] "2068-03-01" "2068-04-01" "2068-05-01" "2068-06-01" "2068-07-01" "2068-08-01" "2068-09-01"
[778] "2068-10-01" "2068-11-01" "2068-12-01" "2069-01-01" "2069-02-01" "2069-03-01" "2069-04-01"
[785] "2069-05-01" "2069-06-01" "2069-07-01" "2069-08-01" "2069-09-01" "2069-10-01" "2069-11-01"
[792] "2069-12-01" "2070-01-01" "2070-02-01" "2070-03-01" "2070-04-01" "2070-05-01" "2070-06-01"
[799] "2070-07-01" "2070-08-01" "2070-09-01" "2070-10-01" "2070-11-01" "2070-12-01" "2071-01-01"
[806] "2071-02-01" "2071-03-01" "2071-04-01" "2071-05-01" "2071-06-01" "2071-07-01" "2071-08-01"
[813] "2071-09-01" "2071-10-01" "2071-11-01" "2071-12-01" "2072-01-01" "2072-02-01" "2072-03-01"
[820] "2072-04-01" "2072-05-01" "2072-06-01" "2072-07-01" "2072-08-01" "2072-09-01" "2072-10-01"
[827] "2072-11-01" "2072-12-01" "2073-01-01" "2073-02-01" "2073-03-01" "2073-04-01" "2073-05-01"
[834] "2073-06-01" "2073-07-01" "2073-08-01" "2073-09-01" "2073-10-01" "2073-11-01" "2073-12-01"
[841] "2074-01-01" "2074-02-01" "2074-03-01" "2074-04-01" "2074-05-01" "2074-06-01" "2074-07-01"
[848] "2074-08-01" "2074-09-01" "2074-10-01" "2074-11-01" "2074-12-01" "2075-01-01" "2075-02-01"
[855] "2075-03-01" "2075-04-01" "2075-05-01" "2075-06-01" "2075-07-01" "2075-08-01" "2075-09-01"
[862] "2075-10-01" "2075-11-01" "2075-12-01" "2076-01-01" "2076-02-01" "2076-03-01" "2076-04-01"
[869] "2076-05-01" "2076-06-01" "2076-07-01" "2076-08-01" "2076-09-01" "2076-10-01" "2076-11-01"
[876] "2076-12-01" "2077-01-01" "2077-02-01" "2077-03-01" "2077-04-01" "2077-05-01" "2077-06-01"
[883] "2077-07-01" "2077-08-01" "2077-09-01" "2077-10-01" "2077-11-01" "2077-12-01" "2078-01-01"
[890] "2078-02-01" "2078-03-01" "2078-04-01" "2078-05-01" "2078-06-01" "2078-07-01" "2078-08-01"
[897] "2078-09-01" "2078-10-01" "2078-11-01" "2078-12-01" "2079-01-01" "2079-02-01" "2079-03-01"
[904] "2079-04-01" "2079-05-01" "2079-06-01" "2079-07-01" "2079-08-01" "2079-09-01" "2079-10-01"
[911] "2079-11-01" "2079-12-01" "2080-01-01" "2080-02-01" "2080-03-01" "2080-04-01" "2080-05-01"
[918] "2080-06-01" "2080-07-01" "2080-08-01" "2080-09-01" "2080-10-01" "2080-11-01" "2080-12-01"
[925] "2081-01-01" "2081-02-01" "2081-03-01" "2081-04-01" "2081-05-01" "2081-06-01" "2081-07-01"
[932] "2081-08-01" "2081-09-01" "2081-10-01" "2081-11-01" "2081-12-01" "2082-01-01" "2082-02-01"
[939] "2082-03-01" "2082-04-01" "2082-05-01" "2082-06-01" "2082-07-01" "2082-08-01" "2082-09-01"
[946] "2082-10-01" "2082-11-01" "2082-12-01" "2083-01-01" "2083-02-01" "2083-03-01" "2083-04-01"
[953] "2083-05-01" "2083-06-01" "2083-07-01" "2083-08-01" "2083-09-01" "2083-10-01" "2083-11-01"
[960] "2083-12-01" "2084-01-01" "2084-02-01" "2084-03-01" "2084-04-01" "2084-05-01" "2084-06-01"
[967] "2084-07-01" "2084-08-01" "2084-09-01" "2084-10-01" "2084-11-01" "2084-12-01" "2085-01-01"
[974] "2085-02-01" "2085-03-01" "2085-04-01" "2085-05-01" "2085-06-01" "2085-07-01" "2085-08-01"
[981] "2085-09-01" "2085-10-01" "2085-11-01" "2085-12-01" "2086-01-01" "2086-02-01" "2086-03-01"
[988] "2086-04-01" "2086-05-01" "2086-06-01" "2086-07-01" "2086-08-01" "2086-09-01" "2086-10-01"
[995] "2086-11-01" "2086-12-01" "2087-01-01" "2087-02-01" "2087-03-01" "2087-04-01"
[ reached getOption("max.print") -- omitted 164 entries ]
time <-as.Date(time)
out <- as.matrix(CRD45_out_response) ## convert from coda to matrix
code from AMW for splitting up this matrix
x_cols <- grep('^x', colnames(out))
obs_cols <- grep('^OBS', colnames(out))
out_x <- out[,x_cols]
out_obs <- out[,obs_cols]
rm(out)
out_x <- t(out_x)
out_x <- as.data.frame(out_x)
out_x2 <- out_x |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*x\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
timesteps <- as.data.frame(time)
timesteps$ind <- as.character(seq(from = 1, to = nrow(timesteps), by = 1))
colnames(timesteps) <- c('time', 'timestep')
ponds <- as.data.frame(colnames(CR_WaterTemp45[-1]))
ponds$ind <- as.character(seq(from = 1, to = nrow(ponds), by = 1))
colnames(ponds) <- c('pond_name', 'pond')
out_x_mapped <- out_x2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pX <- out_x_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
out_obs <- t(out_obs)
out_obs <- as.data.frame(out_obs)
out_obs2 <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
out_obs_mapped <- out_obs2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pOBS <- out_obs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
cowplot::plot_grid(pX, pOBS)
threshold_probs <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(n_iter = dplyr::n(),
n_thresh = sum(val > 20),
prob = (n_thresh / n_iter) * 100)
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
threshold_probs_mapped <- threshold_probs |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond')
threshold_probs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = prob)) +
ggplot2::facet_wrap(~pond_name)
Save these ouputs using a .RDataFile type
save.image(file = "Corr_PooledMod_CRD45.RData")
Load in this file here as long as no code changes are needed above
load(file = "Corr_PooledMod_CRD45.RData")
Creating some nice plots of these
# set up the colors here for the plots below so that each pond has the same color each time
Mod_CRD45 <- out_obs_mapped
str(Mod_CRD45)
gropd_df [10,476 × 6] (S3: grouped_df/tbl_df/tbl/data.frame)
$ pond : chr [1:10476] "1" "1" "1" "1" ...
$ low : Named num [1:10476] -6.549 6.28 -2.714 -0.219 2.658 ...
..- attr(*, "names")= chr [1:10476] "2.5%" "2.5%" "2.5%" "2.5%" ...
$ med : num [1:10476] -0.375 6.28 1.316 3.694 6.286 ...
$ high : Named num [1:10476] 5.95 6.28 5.03 7.62 9.96 ...
..- attr(*, "names")= chr [1:10476] "97.5%" "97.5%" "97.5%" "97.5%" ...
$ time : Date[1:10476], format: "2012-01-01" "2012-10-01" "2012-04-01" ...
$ pond_name: chr [1:10476] "BVS" "BVS" "BVS" "BVS" ...
- attr(*, "groups")= tibble [9 × 2] (S3: tbl_df/tbl/data.frame)
..$ pond : chr [1:9] "1" "2" "3" "4" ...
..$ .rows: list<int> [1:9]
.. ..$ : int [1:1164] 1 2 3 4 5 6 7 8 9 10 ...
.. ..$ : int [1:1164] 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 ...
.. ..$ : int [1:1164] 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 ...
.. ..$ : int [1:1164] 3493 3494 3495 3496 3497 3498 3499 3500 3501 3502 ...
.. ..$ : int [1:1164] 4657 4658 4659 4660 4661 4662 4663 4664 4665 4666 ...
.. ..$ : int [1:1164] 5821 5822 5823 5824 5825 5826 5827 5828 5829 5830 ...
.. ..$ : int [1:1164] 6985 6986 6987 6988 6989 6990 6991 6992 6993 6994 ...
.. ..$ : int [1:1164] 8149 8150 8151 8152 8153 8154 8155 8156 8157 8158 ...
.. ..$ : int [1:1164] 9313 9314 9315 9316 9317 9318 9319 9320 9321 9322 ...
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
write.csv(Mod_CRD45, "Corr_Mod_CRD45.csv")
# Define the custom color palette
CR_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "lightblue")
# Ensure that CR_colors has the same length as the number of unique ponds
unique_ponds <- unique(Mod_CRD45$pond_name)
if (length(CR_colors) != length(unique_ponds)) {
stop("The number of colors in CR_colors does not match the number of unique ponds.")
}
# Create a named vector for CR_colors with pond names
color_mapping <- setNames(CR_colors, unique_ponds)
# Plotting
Ponds_CRD45 <- ggplot(data = Mod_CRD45, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Add the ribbon with colors
geom_line(aes(y = med, color = pond_name)) + # Ensure lines are visible and colored by pond
geom_smooth(aes(y = med, color = "grey"), size = 1) + # Ensure lines are visible and colored by pond
geom_vline(xintercept = as.Date("2020-12-01"), linetype = "dashed", color = "red") + # Add vertical dashed red line
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
xlab("Year") +
ylab("Water Temperature (°C)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Ponds_CRD45
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndPonds_CRD45.png", plot = Ponds_CRD45, width = 8, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# using the same aestetics as above
Trends_CRD45 <- ggplot(data = Mod_CRD45, aes(x = time)) +
geom_smooth(aes(y = med, color = pond_name, fill = pond_name), alpha = 0.15, size = 1) + # Smooth lines with colors by pond_name
scale_color_manual(values = color_mapping) +
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
legend.position = "right",
legend.text = element_text(size = 16), # Change the legend text size here
legend.title = element_text(size = 20)) # Position the legend on the right
Trends_CRD45
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndTrends_CRD45.png", plot = Trends_CRD45, width = 9, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
probs_CRD45 <- threshold_probs_mapped
write.csv(probs_CRD45, "Corr_probs_CRD45.csv")
Probs20_CRD45 <- ggplot(data = probs_CRD45, aes(x = time)) +
geom_line(aes(y = prob, color = pond_name)) + # Ensure lines are visible and colored by pond
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
xlab("Year") +
ylab("Probability") +
xlim(as.Date(c("2020-12-01", "2099-12-01"))) +
ylim(0,50) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 18),
axis.text.x = element_text(angle = 55, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Probs20_CRD45
Warning: Removed 1935 rows containing missing values or values outside the scale range (`geom_line()`).
ggsave("Corr_Probs20_CRD45.png", plot = Probs20_CRD45, width = 9, height = 6)
Warning: Removed 1935 rows containing missing values or values outside the scale range (`geom_line()`).
Mod_CRD45
ModObs <- Mod_CRD45 %>%
filter(time < "2021-01-01")
range(ModObs$time)
[1] "2012-01-01" "2020-12-01"
ModelErrorPlot_CRD45 <- ggplot(data = ModObs, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Adjust fill colors
geom_line(aes(y = med), color = "black") + # Black lines for the median values
facet_wrap(~ pond_name) + # Facet by pond_name
scale_fill_manual(values = color_mapping) + # Apply custom color palette for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
axis.title = element_text(size = 16),
axis.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none" # Angle x-axis labels if needed
)
ModelErrorPlot_CRD45
ggsave("Corr_Error_CRD45.png", plot = ModelErrorPlot_CRD45, width = 9, height = 6)
Clear out the environment and reload the basics before starting the next section
rm(list = ls())
library(tidyverse)
library(rjags)
library(ggplot2)
set.seed(1)
Read in the data here, which is an output of XX code file
Then separate the data into air and water temperature files
# Air and water temperature data: CR loccation, 4.5 scenario
CR_Temps <- read.csv("CR_85SNAPForecast.csv", header=TRUE)
# Air temperature variables
CR_AirTemp85 <- CR_Temps %>%
select(date, Airport, Air_MonthAvg, std, precision, Region)
# Water temperature variables
CR_WaterTemp85 <- CR_Temps %>%
select(date, BVS, CAB, TIN, SQR, CAN, EYS, TIS, WDD, RHM)
# Set up a string of all the dates
time <- CR_WaterTemp85$date
length(time)
[1] 1164
Setting up the random walk model
dlm_pooled <- "
model{
#### Priors
for(p in 1:np){
x[1,p] ~ dnorm(x_ic, tau_ic) # Initial condition of water temperature
}
tau_obs ~ dgamma(a_obs, r_obs) # Prior on observation error
tau_add ~ dgamma(a_add, r_add) # Prior on process error
#### Fixed Effects
beta ~ dmnorm(mu_beta, tau_beta) # Prior on beta coefficients
for(p in 1:np){
int[p] ~ dnorm(mu_int, tau_int) # Prior on pond-specific intercepts
}
#### Data Model
for(t in 1:n){ # loop over all time steps
for(p in 1:np){
OBS[t,p] ~ dnorm(x[t,p], tau_obs) # Observed water temperature is drawn from latent air temperature with observation uncertainty
}
Xf[t] ~ dnorm(muAirTemp[t], tauAirTemp[t]) # Latent air temperature is drawn from mean and precision of forecasated air temperature
}
#### Process Model
for(t in 2:n){ # loop over all time steps except teh first (we defined ic above)
for(p in 1:np){
mu[t,p] <- x[t-1,p] + int[p] + beta[1] * x[t-1,p] + beta[2] * Xf[t] # Mean water temperature is a function of the previous time step and current air temperature
x[t,p] ~ dnorm(mu[t,p], tau_add) # Latent water temperature is drawn from mean water temperature with process uncertainty
}
}
}
"
Define the data and priors for the model
# Empty list
data <- list()
# Water temperature observations
data$OBS <- dplyr::select(CR_WaterTemp85, -date)
# Number of time steps
data$n <- nrow(data$OBS)
# Number of ponds
data$np <- ncol(data$OBS)
# Initial water temperature mean
data$x_ic <- 0.1
# Initial water temperature precision
data$tau_ic = 0.1
# Prior parameters for observation and process uncertainty
data$a_obs = 1
data$r_obs = 1
data$a_add = 1
data$r_add = 1
# Prior parameters for beta coefficients
data$mu_beta <- c(0, 0)
data$tau_beta <- diag(x = c(0.001, 0.001), nrow = 2, ncol = 2)
# Prior parameters for intercept
data$mu_int <- 0
data$tau_int <- 0.001
# Mean air temperature estimate
data$muAirTemp <- CR_AirTemp85$Air_MonthAvg
# Air temperature precision
data$tauAirTemp <- CR_AirTemp85$precision
Create JAGS model with 3 chains
jm <- jags.model(file = textConnection(dlm_pooled), data = data, n.chains = 3)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 448
Unobserved stochastic nodes: 21680
Total graph size: 46573
Initializing model
Posterior samples of parameters
CRD85_out_params <- coda.samples(model = jm,
variable.names = c('beta', 'int',
'tau_add', 'tau_obs'),
n.iter = 50000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Posterior samples of response variables
CRD85_out_response <- coda.samples(model = jm,
variable.names = c('x', 'OBS'),
n.iter = 100000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
A couple quick checks here
# plot(CRD85_out_params)
gelman.diag(CRD85_out_params, confidence = 0.99)
Potential scale reduction factors:
Point est. Upper C.I.
beta[1] 1 1.00
beta[2] 1 1.01
int[1] 1 1.00
int[2] 1 1.02
int[3] 1 1.01
int[4] 1 1.01
int[5] 1 1.01
int[6] 1 1.01
int[7] 1 1.00
int[8] 1 1.01
int[9] 1 1.01
tau_add 1 1.00
tau_obs 1 1.01
Multivariate psrf
1.01
Visualize the output by just looking at the 95% credible interval of the time-series of X’s and compare that to the observed Y’s
Transform the samples back from the log domain to the linear domain
time ## adjust to zoom in and out
[1] "2012-01-01" "2012-02-01" "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01"
[8] "2012-08-01" "2012-09-01" "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01"
[15] "2013-03-01" "2013-04-01" "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01"
[22] "2013-10-01" "2013-11-01" "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01"
[29] "2014-05-01" "2014-06-01" "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01"
[36] "2014-12-01" "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01"
[43] "2015-07-01" "2015-08-01" "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01"
[50] "2016-02-01" "2016-03-01" "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01"
[57] "2016-09-01" "2016-10-01" "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01"
[64] "2017-04-01" "2017-05-01" "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01"
[71] "2017-11-01" "2017-12-01" "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01"
[78] "2018-06-01" "2018-07-01" "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01"
[85] "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01"
[92] "2019-08-01" "2019-09-01" "2019-10-01" "2019-11-01" "2019-12-01" "2012-01-01" "2012-02-01"
[99] "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01" "2012-08-01" "2012-09-01"
[106] "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01" "2013-03-01" "2013-04-01"
[113] "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01" "2013-10-01" "2013-11-01"
[120] "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01" "2014-05-01" "2014-06-01"
[127] "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01" "2014-12-01" "2015-01-01"
[134] "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01" "2015-07-01" "2015-08-01"
[141] "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01" "2016-02-01" "2016-03-01"
[148] "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01"
[155] "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01" "2017-04-01" "2017-05-01"
[162] "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01" "2017-11-01" "2017-12-01"
[169] "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01" "2018-06-01" "2018-07-01"
[176] "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01" "2019-01-01" "2019-02-01"
[183] "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01" "2019-08-01" "2019-09-01"
[190] "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01" "2020-03-01" "2020-04-01"
[197] "2020-05-01" "2020-06-01" "2020-07-01" "2020-08-01" "2020-09-01" "2020-10-01" "2020-11-01"
[204] "2020-12-01" "2021-01-01" "2021-02-01" "2021-03-01" "2021-04-01" "2021-05-01" "2021-06-01"
[211] "2021-07-01" "2021-08-01" "2021-09-01" "2021-10-01" "2021-11-01" "2021-12-01" "2022-01-01"
[218] "2022-02-01" "2022-03-01" "2022-04-01" "2022-05-01" "2022-06-01" "2022-07-01" "2022-08-01"
[225] "2022-09-01" "2022-10-01" "2022-11-01" "2022-12-01" "2023-01-01" "2023-02-01" "2023-03-01"
[232] "2023-04-01" "2023-05-01" "2023-06-01" "2023-07-01" "2023-08-01" "2023-09-01" "2023-10-01"
[239] "2023-11-01" "2023-12-01" "2024-01-01" "2024-02-01" "2024-03-01" "2024-04-01" "2024-05-01"
[246] "2024-06-01" "2024-07-01" "2024-08-01" "2024-09-01" "2024-10-01" "2024-11-01" "2024-12-01"
[253] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01" "2025-06-01" "2025-07-01"
[260] "2025-08-01" "2025-09-01" "2025-10-01" "2025-11-01" "2025-12-01" "2026-01-01" "2026-02-01"
[267] "2026-03-01" "2026-04-01" "2026-05-01" "2026-06-01" "2026-07-01" "2026-08-01" "2026-09-01"
[274] "2026-10-01" "2026-11-01" "2026-12-01" "2027-01-01" "2027-02-01" "2027-03-01" "2027-04-01"
[281] "2027-05-01" "2027-06-01" "2027-07-01" "2027-08-01" "2027-09-01" "2027-10-01" "2027-11-01"
[288] "2027-12-01" "2028-01-01" "2028-02-01" "2028-03-01" "2028-04-01" "2028-05-01" "2028-06-01"
[295] "2028-07-01" "2028-08-01" "2028-09-01" "2028-10-01" "2028-11-01" "2028-12-01" "2029-01-01"
[302] "2029-02-01" "2029-03-01" "2029-04-01" "2029-05-01" "2029-06-01" "2029-07-01" "2029-08-01"
[309] "2029-09-01" "2029-10-01" "2029-11-01" "2029-12-01" "2030-01-01" "2030-02-01" "2030-03-01"
[316] "2030-04-01" "2030-05-01" "2030-06-01" "2030-07-01" "2030-08-01" "2030-09-01" "2030-10-01"
[323] "2030-11-01" "2030-12-01" "2031-01-01" "2031-02-01" "2031-03-01" "2031-04-01" "2031-05-01"
[330] "2031-06-01" "2031-07-01" "2031-08-01" "2031-09-01" "2031-10-01" "2031-11-01" "2031-12-01"
[337] "2032-01-01" "2032-02-01" "2032-03-01" "2032-04-01" "2032-05-01" "2032-06-01" "2032-07-01"
[344] "2032-08-01" "2032-09-01" "2032-10-01" "2032-11-01" "2032-12-01" "2033-01-01" "2033-02-01"
[351] "2033-03-01" "2033-04-01" "2033-05-01" "2033-06-01" "2033-07-01" "2033-08-01" "2033-09-01"
[358] "2033-10-01" "2033-11-01" "2033-12-01" "2034-01-01" "2034-02-01" "2034-03-01" "2034-04-01"
[365] "2034-05-01" "2034-06-01" "2034-07-01" "2034-08-01" "2034-09-01" "2034-10-01" "2034-11-01"
[372] "2034-12-01" "2035-01-01" "2035-02-01" "2035-03-01" "2035-04-01" "2035-05-01" "2035-06-01"
[379] "2035-07-01" "2035-08-01" "2035-09-01" "2035-10-01" "2035-11-01" "2035-12-01" "2036-01-01"
[386] "2036-02-01" "2036-03-01" "2036-04-01" "2036-05-01" "2036-06-01" "2036-07-01" "2036-08-01"
[393] "2036-09-01" "2036-10-01" "2036-11-01" "2036-12-01" "2037-01-01" "2037-02-01" "2037-03-01"
[400] "2037-04-01" "2037-05-01" "2037-06-01" "2037-07-01" "2037-08-01" "2037-09-01" "2037-10-01"
[407] "2037-11-01" "2037-12-01" "2038-01-01" "2038-02-01" "2038-03-01" "2038-04-01" "2038-05-01"
[414] "2038-06-01" "2038-07-01" "2038-08-01" "2038-09-01" "2038-10-01" "2038-11-01" "2038-12-01"
[421] "2039-01-01" "2039-02-01" "2039-03-01" "2039-04-01" "2039-05-01" "2039-06-01" "2039-07-01"
[428] "2039-08-01" "2039-09-01" "2039-10-01" "2039-11-01" "2039-12-01" "2040-01-01" "2040-02-01"
[435] "2040-03-01" "2040-04-01" "2040-05-01" "2040-06-01" "2040-07-01" "2040-08-01" "2040-09-01"
[442] "2040-10-01" "2040-11-01" "2040-12-01" "2041-01-01" "2041-02-01" "2041-03-01" "2041-04-01"
[449] "2041-05-01" "2041-06-01" "2041-07-01" "2041-08-01" "2041-09-01" "2041-10-01" "2041-11-01"
[456] "2041-12-01" "2042-01-01" "2042-02-01" "2042-03-01" "2042-04-01" "2042-05-01" "2042-06-01"
[463] "2042-07-01" "2042-08-01" "2042-09-01" "2042-10-01" "2042-11-01" "2042-12-01" "2043-01-01"
[470] "2043-02-01" "2043-03-01" "2043-04-01" "2043-05-01" "2043-06-01" "2043-07-01" "2043-08-01"
[477] "2043-09-01" "2043-10-01" "2043-11-01" "2043-12-01" "2044-01-01" "2044-02-01" "2044-03-01"
[484] "2044-04-01" "2044-05-01" "2044-06-01" "2044-07-01" "2044-08-01" "2044-09-01" "2044-10-01"
[491] "2044-11-01" "2044-12-01" "2045-01-01" "2045-02-01" "2045-03-01" "2045-04-01" "2045-05-01"
[498] "2045-06-01" "2045-07-01" "2045-08-01" "2045-09-01" "2045-10-01" "2045-11-01" "2045-12-01"
[505] "2046-01-01" "2046-02-01" "2046-03-01" "2046-04-01" "2046-05-01" "2046-06-01" "2046-07-01"
[512] "2046-08-01" "2046-09-01" "2046-10-01" "2046-11-01" "2046-12-01" "2047-01-01" "2047-02-01"
[519] "2047-03-01" "2047-04-01" "2047-05-01" "2047-06-01" "2047-07-01" "2047-08-01" "2047-09-01"
[526] "2047-10-01" "2047-11-01" "2047-12-01" "2048-01-01" "2048-02-01" "2048-03-01" "2048-04-01"
[533] "2048-05-01" "2048-06-01" "2048-07-01" "2048-08-01" "2048-09-01" "2048-10-01" "2048-11-01"
[540] "2048-12-01" "2049-01-01" "2049-02-01" "2049-03-01" "2049-04-01" "2049-05-01" "2049-06-01"
[547] "2049-07-01" "2049-08-01" "2049-09-01" "2049-10-01" "2049-11-01" "2049-12-01" "2050-01-01"
[554] "2050-02-01" "2050-03-01" "2050-04-01" "2050-05-01" "2050-06-01" "2050-07-01" "2050-08-01"
[561] "2050-09-01" "2050-10-01" "2050-11-01" "2050-12-01" "2051-01-01" "2051-02-01" "2051-03-01"
[568] "2051-04-01" "2051-05-01" "2051-06-01" "2051-07-01" "2051-08-01" "2051-09-01" "2051-10-01"
[575] "2051-11-01" "2051-12-01" "2052-01-01" "2052-02-01" "2052-03-01" "2052-04-01" "2052-05-01"
[582] "2052-06-01" "2052-07-01" "2052-08-01" "2052-09-01" "2052-10-01" "2052-11-01" "2052-12-01"
[589] "2053-01-01" "2053-02-01" "2053-03-01" "2053-04-01" "2053-05-01" "2053-06-01" "2053-07-01"
[596] "2053-08-01" "2053-09-01" "2053-10-01" "2053-11-01" "2053-12-01" "2054-01-01" "2054-02-01"
[603] "2054-03-01" "2054-04-01" "2054-05-01" "2054-06-01" "2054-07-01" "2054-08-01" "2054-09-01"
[610] "2054-10-01" "2054-11-01" "2054-12-01" "2055-01-01" "2055-02-01" "2055-03-01" "2055-04-01"
[617] "2055-05-01" "2055-06-01" "2055-07-01" "2055-08-01" "2055-09-01" "2055-10-01" "2055-11-01"
[624] "2055-12-01" "2056-01-01" "2056-02-01" "2056-03-01" "2056-04-01" "2056-05-01" "2056-06-01"
[631] "2056-07-01" "2056-08-01" "2056-09-01" "2056-10-01" "2056-11-01" "2056-12-01" "2057-01-01"
[638] "2057-02-01" "2057-03-01" "2057-04-01" "2057-05-01" "2057-06-01" "2057-07-01" "2057-08-01"
[645] "2057-09-01" "2057-10-01" "2057-11-01" "2057-12-01" "2058-01-01" "2058-02-01" "2058-03-01"
[652] "2058-04-01" "2058-05-01" "2058-06-01" "2058-07-01" "2058-08-01" "2058-09-01" "2058-10-01"
[659] "2058-11-01" "2058-12-01" "2059-01-01" "2059-02-01" "2059-03-01" "2059-04-01" "2059-05-01"
[666] "2059-06-01" "2059-07-01" "2059-08-01" "2059-09-01" "2059-10-01" "2059-11-01" "2059-12-01"
[673] "2060-01-01" "2060-02-01" "2060-03-01" "2060-04-01" "2060-05-01" "2060-06-01" "2060-07-01"
[680] "2060-08-01" "2060-09-01" "2060-10-01" "2060-11-01" "2060-12-01" "2061-01-01" "2061-02-01"
[687] "2061-03-01" "2061-04-01" "2061-05-01" "2061-06-01" "2061-07-01" "2061-08-01" "2061-09-01"
[694] "2061-10-01" "2061-11-01" "2061-12-01" "2062-01-01" "2062-02-01" "2062-03-01" "2062-04-01"
[701] "2062-05-01" "2062-06-01" "2062-07-01" "2062-08-01" "2062-09-01" "2062-10-01" "2062-11-01"
[708] "2062-12-01" "2063-01-01" "2063-02-01" "2063-03-01" "2063-04-01" "2063-05-01" "2063-06-01"
[715] "2063-07-01" "2063-08-01" "2063-09-01" "2063-10-01" "2063-11-01" "2063-12-01" "2064-01-01"
[722] "2064-02-01" "2064-03-01" "2064-04-01" "2064-05-01" "2064-06-01" "2064-07-01" "2064-08-01"
[729] "2064-09-01" "2064-10-01" "2064-11-01" "2064-12-01" "2065-01-01" "2065-02-01" "2065-03-01"
[736] "2065-04-01" "2065-05-01" "2065-06-01" "2065-07-01" "2065-08-01" "2065-09-01" "2065-10-01"
[743] "2065-11-01" "2065-12-01" "2066-01-01" "2066-02-01" "2066-03-01" "2066-04-01" "2066-05-01"
[750] "2066-06-01" "2066-07-01" "2066-08-01" "2066-09-01" "2066-10-01" "2066-11-01" "2066-12-01"
[757] "2067-01-01" "2067-02-01" "2067-03-01" "2067-04-01" "2067-05-01" "2067-06-01" "2067-07-01"
[764] "2067-08-01" "2067-09-01" "2067-10-01" "2067-11-01" "2067-12-01" "2068-01-01" "2068-02-01"
[771] "2068-03-01" "2068-04-01" "2068-05-01" "2068-06-01" "2068-07-01" "2068-08-01" "2068-09-01"
[778] YF_WaterTemp85 "2068-11-01" "2068-12-01" "2069-01-01" "2069-02-01" "2069-03-01" "2069-04-01"
[785] "2069-05-01" "2069-06-01" "2069-07-01" "2069-08-01" "2069-09-01" "2069-10-01" "2069-11-01"
[792] "2069-12-01" "2070-01-01" "2070-02-01" "2070-03-01" "2070-04-01" "2070-05-01" "2070-06-01"
[799] "2070-07-01" "2070-08-01" "2070-09-01" "2070-10-01" "2070-11-01" "2070-12-01" "2071-01-01"
[806] "2071-02-01" "2071-03-01" "2071-04-01" "2071-05-01" "2071-06-01" "2071-07-01" "2071-08-01"
[813] "2071-09-01" "2071-10-01" "2071-11-01" "2071-12-01" "2072-01-01" "2072-02-01" "2072-03-01"
[820] "2072-04-01" "2072-05-01" "2072-06-01" "2072-07-01" "2072-08-01" "2072-09-01" "2072-10-01"
[827] "2072-11-01" "2072-12-01" "2073-01-01" "2073-02-01" "2073-03-01" "2073-04-01" "2073-05-01"
[834] "2073-06-01" "2073-07-01" "2073-08-01" "2073-09-01" "2073-10-01" "2073-11-01" "2073-12-01"
[841] "2074-01-01" "2074-02-01" "2074-03-01" "2074-04-01" "2074-05-01" "2074-06-01" "2074-07-01"
[848] "2074-08-01" "2074-09-01" "2074-10-01" "2074-11-01" "2074-12-01" "2075-01-01" "2075-02-01"
[855] "2075-03-01" "2075-04-01" "2075-05-01" "2075-06-01" "2075-07-01" "2075-08-01" "2075-09-01"
[862] "2075-10-01" "2075-11-01" "2075-12-01" "2076-01-01" "2076-02-01" "2076-03-01" "2076-04-01"
[869] "2076-05-01" "2076-06-01" "2076-07-01" "2076-08-01" "2076-09-01" "2076-10-01" "2076-11-01"
[876] "2076-12-01" "2077-01-01" "2077-02-01" "2077-03-01" "2077-04-01" "2077-05-01" "2077-06-01"
[883] "2077-07-01" "2077-08-01" "2077-09-01" "2077-10-01" "2077-11-01" "2077-12-01" "2078-01-01"
[890] "2078-02-01" "2078-03-01" "2078-04-01" "2078-05-01" "2078-06-01" "2078-07-01" "2078-08-01"
[897] "2078-09-01" "2078-10-01" "2078-11-01" "2078-12-01" "2079-01-01" "2079-02-01" "2079-03-01"
[904] "2079-04-01" "2079-05-01" "2079-06-01" "2079-07-01" "2079-08-01" "2079-09-01" "2079-10-01"
[911] "2079-11-01" "2079-12-01" "2080-01-01" "2080-02-01" "2080-03-01" "2080-04-01" "2080-05-01"
[918] "2080-06-01" "2080-07-01" "2080-08-01" "2080-09-01" "2080-10-01" "2080-11-01" "2080-12-01"
[925] "2081-01-01" "2081-02-01" "2081-03-01" "2081-04-01" "2081-05-01" "2081-06-01" "2081-07-01"
[932] "2081-08-01" "2081-09-01" "2081-10-01" "2081-11-01" "2081-12-01" "2082-01-01" "2082-02-01"
[939] "2082-03-01" "2082-04-01" "2082-05-01" "2082-06-01" "2082-07-01" "2082-08-01" "2082-09-01"
[946] "2082-10-01" "2082-11-01" "2082-12-01" "2083-01-01" "2083-02-01" "2083-03-01" "2083-04-01"
[953] "2083-05-01" "2083-06-01" "2083-07-01" "2083-08-01" "2083-09-01" "2083-10-01" "2083-11-01"
[960] "2083-12-01" "2084-01-01" "2084-02-01" "2084-03-01" "2084-04-01" "2084-05-01" "2084-06-01"
[967] "2084-07-01" "2084-08-01" "2084-09-01" "2084-10-01" "2084-11-01" "2084-12-01" "2085-01-01"
[974] "2085-02-01" "2085-03-01" "2085-04-01" "2085-05-01" "2085-06-01" "2085-07-01" "2085-08-01"
[981] "2085-09-01" "2085-10-01" "2085-11-01" "2085-12-01" "2086-01-01" "2086-02-01" "2086-03-01"
[988] "2086-04-01" "2086-05-01" "2086-06-01" "2086-07-01" "2086-08-01" "2086-09-01" "2086-10-01"
[995] "2086-11-01" "2086-12-01" "2087-01-01" "2087-02-01" "2087-03-01" "2087-04-01"
[ reached getOption("max.print") -- omitted 164 entries ]
time <-as.Date(time)
out <- as.matrix(CRD85_out_response) ## convert from coda to matrix
code from AMW for splitting up this matrix
x_cols <- grep('^x', colnames(out))
obs_cols <- grep('^OBS', colnames(out))
out_x <- out[,x_cols]
out_obs <- out[,obs_cols]
rm(out)
out_x <- t(out_x)
out_x <- as.data.frame(out_x)
out_x2 <- out_x |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*x\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
timesteps <- as.data.frame(time)
timesteps$ind <- as.character(seq(from = 1, to = nrow(timesteps), by = 1))
colnames(timesteps) <- c('time', 'timestep')
ponds <- as.data.frame(colnames(CR_WaterTemp85[-1]))
ponds$ind <- as.character(seq(from = 1, to = nrow(ponds), by = 1))
colnames(ponds) <- c('pond_name', 'pond')
out_x_mapped <- out_x2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pX <- out_x_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
out_obs <- t(out_obs)
out_obs <- as.data.frame(out_obs)
out_obs2 <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
out_obs_mapped <- out_obs2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pOBS <- out_obs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
cowplot::plot_grid(pX, pOBS)
threshold_probs <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(n_iter = dplyr::n(),
n_thresh = sum(val > 20),
prob = (n_thresh / n_iter) * 100)
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
threshold_probs_mapped <- threshold_probs |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond')
threshold_probs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = prob)) +
ggplot2::facet_wrap(~pond_name)
Save these ouputs using a .RDataFile type
save.image(file = "Corr_PooledMod_CRD85.RData")
#Currently this has all the information in it. Eventually could limit this to be just the relavent information for the remaining code below
Load in this file here as long as no code changes are needed above
load(file = "Corr_PooledMod_CRD85.RData")
Creating some nice plots of these
# set up the colors here for the plots below so that each pond has the same color each time
Mod_CRD85 <- out_obs_mapped
str(Mod_CRD85)
gropd_df [10,476 × 6] (S3: grouped_df/tbl_df/tbl/data.frame)
$ pond : chr [1:10476] "1" "1" "1" "1" ...
$ low : Named num [1:10476] -6.53 6.28 -3.61 1.54 3.74 ...
..- attr(*, "names")= chr [1:10476] "2.5%" "2.5%" "2.5%" "2.5%" ...
$ med : num [1:10476] -0.351 6.28 0.736 5.619 8.518 ...
$ high : Named num [1:10476] 5.77 6.28 4.8 9.64 13.36 ...
..- attr(*, "names")= chr [1:10476] "97.5%" "97.5%" "97.5%" "97.5%" ...
$ time : Date[1:10476], format: "2012-01-01" "2012-10-01" "2012-04-01" ...
$ pond_name: chr [1:10476] "BVS" "BVS" "BVS" "BVS" ...
- attr(*, "groups")= tibble [9 × 2] (S3: tbl_df/tbl/data.frame)
..$ pond : chr [1:9] "1" "2" "3" "4" ...
..$ .rows: list<int> [1:9]
.. ..$ : int [1:1164] 1 2 3 4 5 6 7 8 9 10 ...
.. ..$ : int [1:1164] 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 ...
.. ..$ : int [1:1164] 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 ...
.. ..$ : int [1:1164] 3493 3494 3495 3496 3497 3498 3499 3500 3501 3502 ...
.. ..$ : int [1:1164] 4657 4658 4659 4660 4661 4662 4663 4664 4665 4666 ...
.. ..$ : int [1:1164] 5821 5822 5823 5824 5825 5826 5827 5828 5829 5830 ...
.. ..$ : int [1:1164] 6985 6986 6987 6988 6989 6990 6991 6992 6993 6994 ...
.. ..$ : int [1:1164] 8149 8150 8151 8152 8153 8154 8155 8156 8157 8158 ...
.. ..$ : int [1:1164] 9313 9314 9315 9316 9317 9318 9319 9320 9321 9322 ...
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
write.csv(Mod_CRD85, "Corr_Mod_CRD85.csv")
# Define the custom color palette
CR_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "lightblue")
# Ensure that CR_colors has the same length as the number of unique ponds
unique_ponds <- unique(Mod_CRD85$pond_name)
if (length(CR_colors) != length(unique_ponds)) {
stop("The number of colors in CR_colors does not match the number of unique ponds.")
}
# Create a named vector for CR_colors with pond names
color_mapping <- setNames(CR_colors, unique_ponds)
# Plotting
Ponds_CRD85 <- ggplot(data = Mod_CRD85, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Add the ribbon with colors
geom_line(aes(y = med, color = pond_name)) + # Ensure lines are visible and colored by pond
geom_smooth(aes(y = med, color = "grey"), size = 1) + # Ensure lines are visible and colored by pond
geom_vline(xintercept = as.Date("2020-12-01"), linetype = "dashed", color = "red") + # Add vertical dashed red line
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
xlab("Year") +
ylab("Water Temperature (°C)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Ponds_CRD85
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndPonds_CRD85.png", plot = Ponds_CRD85, width = 8, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# using the same aestetics as above
Trends_CRD85 <- ggplot(data = Mod_CRD85, aes(x = time)) +
geom_smooth(aes(y = med, color = pond_name, fill = pond_name), alpha = 0.15, size = 1) + # Smooth lines with colors by pond_name
scale_color_manual(values = color_mapping) +
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
legend.position = "right",
legend.text = element_text(size = 16), # Change the legend text size here
legend.title = element_text(size = 20))
Trends_CRD85
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndTrends_CRD85.png", plot = Trends_CRD85, width = 9, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
probs_CRD85 <- threshold_probs_mapped
write.csv(probs_CRD85, "probs_CRD85.csv")
Probs20_CRD85 <- ggplot(data = probs_CRD85, aes(x = time)) +
geom_line(aes(y = prob, color = pond_name)) + # Ensure lines are visible and colored by pond
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
xlab("Year") +
ylab("Probability") +
xlim(as.Date(c("2020-12-01", "2099-12-01"))) +
ylim(0,50) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 18),
axis.text.x = element_text(angle = 55, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Probs20_CRD85
Warning: Removed 1935 rows containing missing values or values outside the scale range (`geom_line()`).
ggsave("Corr_Probs20_CRD85.png", plot = Probs20_CRD85, width = 9, height = 6)
Warning: Removed 1935 rows containing missing values or values outside the scale range (`geom_line()`).
Mod_CRD85
ModObs <- Mod_CRD85 %>%
filter(time < "2021-01-01")
range(ModObs$time)
[1] "2012-01-01" "2020-12-01"
ModelErrorPlot_CRD85 <- ggplot(data = ModObs, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Adjust fill colors
geom_line(aes(y = med), color = "black") + # Black lines for the median values
facet_wrap(~ pond_name) + # Facet by pond_name
scale_fill_manual(values = color_mapping) + # Apply custom color palette for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
ModelErrorPlot_CRD85
ggsave("Corr_Error_CRD85.png", plot = ModelErrorPlot_CRD85, width = 9, height = 6)
Mod_CRD85
# Filter this to only the observational data
obs_CRD85 <- Mod_CRD85 %>%
filter(time < as.Date("2020-01-01")) %>%
ungroup() %>%
as.data.frame()
# Convert the 'date' column in CR_WaterTemp85 to Date type and rename it to 'time'
CR_WaterTemp85 <- CR_WaterTemp85 %>%
mutate(time = as.Date(date))
# BVS ####
BVS_CRD <- obs_CRD85 %>%
filter(pond_name == "BVS")
# Perform the merge to combine the data frames based on the time/date columns
BVS_combined <- BVS_CRD %>%
select(time, med) %>%
left_join(CR_WaterTemp85, by = "time")
Warning in left_join(., CR_WaterTemp85, by = "time") :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 10 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this
warning.
BVS_NSE <- NSE(BVS_combined$med, BVS_combined$BVS, na.rm = TRUE)
BVS_NSE
[1] 0.7038907
# CAB ####
CAB_CRD <- obs_CRD85 %>%
filter(pond_name == "CAB")
# Perform the merge to combine the data frames based on the time/date columns
CAB_combined <- CAB_CRD %>%
select(time, med) %>%
left_join(CR_WaterTemp85, by = "time")
Warning in left_join(., CR_WaterTemp85, by = "time") :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 10 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this
warning.
CAB_NSE <- NSE(CAB_combined$med, CAB_combined$CAB, na.rm = TRUE)
CAB_NSE
[1] 0.8513042
# CAN ####
CAN_CRD <- obs_CRD85 %>%
filter(pond_name == "CAN")
# Perform the merge to combine the data frames based on the time/date columns
CAN_combined <- CAN_CRD %>%
select(time, med) %>%
left_join(CR_WaterTemp85, by = "time")
Warning in left_join(., CR_WaterTemp85, by = "time") :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 10 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this
warning.
CAN_NSE <- NSE(CAN_combined$med, CAN_combined$CAN, na.rm = TRUE)
CAN_NSE
[1] 0.8422728
# CAN ####
CAN_CRD <- obs_CRD85 %>%
filter(pond_name == "CAN")
# Perform the merge to combine the data frames based on the time/date columns
CAN_combined <- CAN_CRD %>%
select(time, med) %>%
left_join(CR_WaterTemp85, by = "time")
Warning in left_join(., CR_WaterTemp85, by = "time") :
Detected an unexpected many-to-many relationship between `x` and `y`.
ℹ Row 1 of `x` matches multiple rows in `y`.
ℹ Row 10 of `y` matches multiple rows in `x`.
ℹ If a many-to-many relationship is expected, set `relationship = "many-to-many"` to silence this
warning.
CAN_NSE <- NSE(CAN_combined$med, CAN_combined$CAN, na.rm = TRUE)
CAN_NSE
[1] 0.8422728
Clear out the environment and reload the basics before starting the next section
rm(list = ls())
library(tidyverse)
library(rjags)
library(ggplot2)
set.seed(1)
Read in the data here, which is an output of XX code file
Then separate the data into air and water temperature files
# Air and water temperature data: YF loccation, 4.5 scenario
YF_Temps <- read.csv("BiasCorrection/Geo_BiasCorrection/DataFiles/SNAP_YF45_corrected_monthly.csv", header=TRUE)
# Air temperature variables
YF_AirTemp45 <- YF_Temps %>%
select(date, Airport, Air_MonthAvg, std, precision, Region)
Error in `select()`:
! Can't select columns that don't exist.
✖ Column `Airport` doesn't exist.
Backtrace:
1. YF_Temps %>% ...
3. dplyr:::select.data.frame(., date, Airport, Air_MonthAvg, std, precision, Region)
Setting up the random walk model
dlm_pooled <- "
model{
#### Priors
for(p in 1:np){
x[1,p] ~ dnorm(x_ic, tau_ic) # Initial condition of water temperature
}
tau_obs ~ dgamma(a_obs, r_obs) # Prior on observation error
tau_add ~ dgamma(a_add, r_add) # Prior on process error
#### Fixed Effects
beta ~ dmnorm(mu_beta, tau_beta) # Prior on beta coefficients
for(p in 1:np){
int[p] ~ dnorm(mu_int, tau_int) # Prior on pond-specific intercepts
}
#### Data Model
for(t in 1:n){ # loop over all time steps
for(p in 1:np){
OBS[t,p] ~ dnorm(x[t,p], tau_obs) # Observed water temperature is drawn from latent air temperature with observation uncertainty
}
Xf[t] ~ dnorm(muAirTemp[t], tauAirTemp[t]) # Latent air temperature is drawn from mean and precision of forecasated air temperature
}
#### Process Model
for(t in 2:n){ # loop over all time steps except teh first (we defined ic above)
for(p in 1:np){
mu[t,p] <- x[t-1,p] + int[p] + beta[1] * x[t-1,p] + beta[2] * Xf[t] # Mean water temperature is a function of the previous time step and current air temperature
x[t,p] ~ dnorm(mu[t,p], tau_add) # Latent water temperature is drawn from mean water temperature with process uncertainty
}
}
}
"
Define the data and priors for the model
# Empty list
data <- list()
# Water temperature observations
data$OBS <- dplyr::select(YF_WaterTemp45, -date)
# Number of time steps
data$n <- nrow(data$OBS)
# Number of ponds
data$np <- ncol(data$OBS)
# Initial water temperature mean
data$x_ic <- 0.1
# Initial water temperature precision
data$tau_ic = 0.1
# Prior parameters for observation and process uncertainty
data$a_obs = 1
data$r_obs = 1
data$a_add = 1
data$r_add = 1
# Prior parameters for beta coefficients
data$mu_beta <- c(0, 0)
data$tau_beta <- diag(x = c(0.001, 0.001), nrow = 2, ncol = 2)
# Prior parameters for intercept
data$mu_int <- 0
data$tau_int <- 0.001
# Mean air temperature estimate
data$muAirTemp <- YF_AirTemp45$Air_MonthAvg
# Air temperature precision
data$tauAirTemp <- YF_AirTemp45$precision
Create JAGS model with 3 chains
jm <- jags.model(file = textConnection(dlm_pooled), data = data, n.chains = 3)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 456
Unobserved stochastic nodes: 26330
Total graph size: 55883
Initializing model
Posterior samples of parameters
YF45_out_params <- coda.samples(model = jm,
variable.names = c('beta', 'int',
'tau_add', 'tau_obs'),
n.iter = 100000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Posterior samples of response variables
YF45_out_response <- coda.samples(model = jm,
variable.names = c('x', 'OBS'),
n.iter = 100000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
A couple quick checks here
# plot(YF45_out_params)
gelman.diag(YF45_out_params, confidence = 0.99)
Potential scale reduction factors:
Point est. Upper C.I.
beta[1] 1.00 1.02
beta[2] 1.00 1.01
int[1] 1.01 1.03
int[2] 1.01 1.02
int[3] 1.01 1.03
int[4] 1.01 1.03
int[5] 1.01 1.03
int[6] 1.01 1.03
int[7] 1.00 1.02
int[8] 1.01 1.02
int[9] 1.01 1.03
int[10] 1.01 1.03
int[11] 1.01 1.02
tau_add 1.00 1.01
tau_obs 1.00 1.00
Multivariate psrf
1.01
Visualize the output by just looking at the 95% credible interval of the time-series of X’s and compare that to the observed Y’s
Transform the samples back from the log domain to the linear domain
time ## adjust to zoom in and out
[1] "2012-01-01" "2012-02-01" "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01"
[8] "2012-08-01" "2012-09-01" "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01"
[15] "2013-03-01" "2013-04-01" "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01"
[22] "2013-10-01" "2013-11-01" "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01"
[29] "2014-05-01" "2014-06-01" "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01"
[36] "2014-12-01" "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01"
[43] "2015-07-01" "2015-08-01" "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01"
[50] "2016-02-01" "2016-03-01" "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01"
[57] "2016-09-01" "2016-10-01" "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01"
[64] "2017-04-01" "2017-05-01" "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01"
[71] "2017-11-01" "2017-12-01" "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01"
[78] "2018-06-01" "2018-07-01" "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01"
[85] "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01"
[92] "2019-08-01" "2019-09-01" "2019-10-01" "2019-11-01" "2019-12-01" "2012-01-01" "2012-02-01"
[99] "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01" "2012-08-01" "2012-09-01"
[106] "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01" "2013-03-01" "2013-04-01"
[113] "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01" "2013-10-01" "2013-11-01"
[120] "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01" "2014-05-01" "2014-06-01"
[127] "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01" "2014-12-01" "2015-01-01"
[134] "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01" "2015-07-01" "2015-08-01"
[141] "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01" "2016-02-01" "2016-03-01"
[148] "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01"
[155] "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01" "2017-04-01" "2017-05-01"
[162] "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01" "2017-11-01" "2017-12-01"
[169] "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01" "2018-06-01" "2018-07-01"
[176] "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01" "2019-01-01" "2019-02-01"
[183] "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01" "2019-08-01" "2019-09-01"
[190] "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01" "2020-03-01" "2020-04-01"
[197] "2020-05-01" "2020-06-01" "2020-07-01" "2020-08-01" "2020-09-01" "2020-10-01" "2020-11-01"
[204] "2020-12-01" "2021-01-01" "2021-02-01" "2021-03-01" "2021-04-01" "2021-05-01" "2021-06-01"
[211] "2021-07-01" "2021-08-01" "2021-09-01" "2021-10-01" "2021-11-01" "2021-12-01" "2022-01-01"
[218] "2022-02-01" "2022-03-01" "2022-04-01" "2022-05-01" "2022-06-01" "2022-07-01" "2022-08-01"
[225] "2022-09-01" "2022-10-01" "2022-11-01" "2022-12-01" "2023-01-01" "2023-02-01" "2023-03-01"
[232] "2023-04-01" "2023-05-01" "2023-06-01" "2023-07-01" "2023-08-01" "2023-09-01" "2023-10-01"
[239] "2023-11-01" "2023-12-01" "2024-01-01" "2024-02-01" "2024-03-01" "2024-04-01" "2024-05-01"
[246] "2024-06-01" "2024-07-01" "2024-08-01" "2024-09-01" "2024-10-01" "2024-11-01" "2024-12-01"
[253] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01" "2025-06-01" "2025-07-01"
[260] "2025-08-01" "2025-09-01" "2025-10-01" "2025-11-01" "2025-12-01" "2026-01-01" "2026-02-01"
[267] "2026-03-01" "2026-04-01" "2026-05-01" "2026-06-01" "2026-07-01" "2026-08-01" "2026-09-01"
[274] "2026-10-01" "2026-11-01" "2026-12-01" "2027-01-01" "2027-02-01" "2027-03-01" "2027-04-01"
[281] "2027-05-01" "2027-06-01" "2027-07-01" "2027-08-01" "2027-09-01" "2027-10-01" "2027-11-01"
[288] "2027-12-01" "2028-01-01" "2028-02-01" "2028-03-01" "2028-04-01" "2028-05-01" "2028-06-01"
[295] "2028-07-01" "2028-08-01" "2028-09-01" "2028-10-01" "2028-11-01" "2028-12-01" "2029-01-01"
[302] "2029-02-01" "2029-03-01" "2029-04-01" "2029-05-01" "2029-06-01" "2029-07-01" "2029-08-01"
[309] "2029-09-01" "2029-10-01" "2029-11-01" "2029-12-01" "2030-01-01" "2030-02-01" "2030-03-01"
[316] "2030-04-01" "2030-05-01" "2030-06-01" "2030-07-01" "2030-08-01" "2030-09-01" "2030-10-01"
[323] "2030-11-01" "2030-12-01" "2031-01-01" "2031-02-01" "2031-03-01" "2031-04-01" "2031-05-01"
[330] "2031-06-01" "2031-07-01" "2031-08-01" "2031-09-01" "2031-10-01" "2031-11-01" "2031-12-01"
[337] "2032-01-01" "2032-02-01" "2032-03-01" "2032-04-01" "2032-05-01" "2032-06-01" "2032-07-01"
[344] "2032-08-01" "2032-09-01" "2032-10-01" "2032-11-01" "2032-12-01" "2033-01-01" "2033-02-01"
[351] "2033-03-01" "2033-04-01" "2033-05-01" "2033-06-01" "2033-07-01" "2033-08-01" "2033-09-01"
[358] "2033-10-01" "2033-11-01" "2033-12-01" "2034-01-01" "2034-02-01" "2034-03-01" "2034-04-01"
[365] "2034-05-01" "2034-06-01" "2034-07-01" "2034-08-01" "2034-09-01" "2034-10-01" "2034-11-01"
[372] "2034-12-01" "2035-01-01" "2035-02-01" "2035-03-01" "2035-04-01" "2035-05-01" "2035-06-01"
[379] "2035-07-01" "2035-08-01" "2035-09-01" "2035-10-01" "2035-11-01" "2035-12-01" "2036-01-01"
[386] "2036-02-01" "2036-03-01" "2036-04-01" "2036-05-01" "2036-06-01" "2036-07-01" "2036-08-01"
[393] "2036-09-01" "2036-10-01" "2036-11-01" "2036-12-01" "2037-01-01" "2037-02-01" "2037-03-01"
[400] "2037-04-01" "2037-05-01" "2037-06-01" "2037-07-01" "2037-08-01" "2037-09-01" "2037-10-01"
[407] "2037-11-01" "2037-12-01" "2038-01-01" "2038-02-01" "2038-03-01" "2038-04-01" "2038-05-01"
[414] "2038-06-01" "2038-07-01" "2038-08-01" "2038-09-01" "2038-10-01" "2038-11-01" "2038-12-01"
[421] "2039-01-01" "2039-02-01" "2039-03-01" "2039-04-01" "2039-05-01" "2039-06-01" "2039-07-01"
[428] "2039-08-01" "2039-09-01" "2039-10-01" "2039-11-01" "2039-12-01" "2040-01-01" "2040-02-01"
[435] "2040-03-01" "2040-04-01" "2040-05-01" "2040-06-01" "2040-07-01" "2040-08-01" "2040-09-01"
[442] "2040-10-01" "2040-11-01" "2040-12-01" "2041-01-01" "2041-02-01" "2041-03-01" "2041-04-01"
[449] "2041-05-01" "2041-06-01" "2041-07-01" "2041-08-01" "2041-09-01" "2041-10-01" "2041-11-01"
[456] "2041-12-01" "2042-01-01" "2042-02-01" "2042-03-01" "2042-04-01" "2042-05-01" "2042-06-01"
[463] "2042-07-01" "2042-08-01" "2042-09-01" "2042-10-01" "2042-11-01" "2042-12-01" "2043-01-01"
[470] "2043-02-01" "2043-03-01" "2043-04-01" "2043-05-01" "2043-06-01" "2043-07-01" "2043-08-01"
[477] "2043-09-01" "2043-10-01" "2043-11-01" "2043-12-01" "2044-01-01" "2044-02-01" "2044-03-01"
[484] "2044-04-01" "2044-05-01" "2044-06-01" "2044-07-01" "2044-08-01" "2044-09-01" "2044-10-01"
[491] "2044-11-01" "2044-12-01" "2045-01-01" "2045-02-01" "2045-03-01" "2045-04-01" "2045-05-01"
[498] "2045-06-01" "2045-07-01" "2045-08-01" "2045-09-01" "2045-10-01" "2045-11-01" "2045-12-01"
[505] "2046-01-01" "2046-02-01" "2046-03-01" "2046-04-01" "2046-05-01" "2046-06-01" "2046-07-01"
[512] "2046-08-01" "2046-09-01" "2046-10-01" "2046-11-01" "2046-12-01" "2047-01-01" "2047-02-01"
[519] "2047-03-01" "2047-04-01" "2047-05-01" "2047-06-01" "2047-07-01" "2047-08-01" "2047-09-01"
[526] "2047-10-01" "2047-11-01" "2047-12-01" "2048-01-01" "2048-02-01" "2048-03-01" "2048-04-01"
[533] "2048-05-01" "2048-06-01" "2048-07-01" "2048-08-01" "2048-09-01" "2048-10-01" "2048-11-01"
[540] "2048-12-01" "2049-01-01" "2049-02-01" "2049-03-01" "2049-04-01" "2049-05-01" "2049-06-01"
[547] "2049-07-01" "2049-08-01" "2049-09-01" "2049-10-01" "2049-11-01" "2049-12-01" "2050-01-01"
[554] "2050-02-01" "2050-03-01" "2050-04-01" "2050-05-01" "2050-06-01" "2050-07-01" "2050-08-01"
[561] "2050-09-01" "2050-10-01" "2050-11-01" "2050-12-01" "2051-01-01" "2051-02-01" "2051-03-01"
[568] "2051-04-01" "2051-05-01" "2051-06-01" "2051-07-01" "2051-08-01" "2051-09-01" "2051-10-01"
[575] "2051-11-01" "2051-12-01" "2052-01-01" "2052-02-01" "2052-03-01" "2052-04-01" "2052-05-01"
[582] "2052-06-01" "2052-07-01" "2052-08-01" "2052-09-01" "2052-10-01" "2052-11-01" "2052-12-01"
[589] "2053-01-01" "2053-02-01" "2053-03-01" "2053-04-01" "2053-05-01" "2053-06-01" "2053-07-01"
[596] "2053-08-01" "2053-09-01" "2053-10-01" "2053-11-01" "2053-12-01" "2054-01-01" "2054-02-01"
[603] "2054-03-01" "2054-04-01" "2054-05-01" "2054-06-01" "2054-07-01" "2054-08-01" "2054-09-01"
[610] "2054-10-01" "2054-11-01" "2054-12-01" "2055-01-01" "2055-02-01" "2055-03-01" "2055-04-01"
[617] "2055-05-01" "2055-06-01" "2055-07-01" "2055-08-01" "2055-09-01" "2055-10-01" "2055-11-01"
[624] "2055-12-01" "2056-01-01" "2056-02-01" "2056-03-01" "2056-04-01" "2056-05-01" "2056-06-01"
[631] "2056-07-01" "2056-08-01" "2056-09-01" "2056-10-01" "2056-11-01" "2056-12-01" "2057-01-01"
[638] "2057-02-01" "2057-03-01" "2057-04-01" "2057-05-01" "2057-06-01" "2057-07-01" "2057-08-01"
[645] "2057-09-01" "2057-10-01" "2057-11-01" "2057-12-01" "2058-01-01" "2058-02-01" "2058-03-01"
[652] "2058-04-01" "2058-05-01" "2058-06-01" "2058-07-01" "2058-08-01" "2058-09-01" "2058-10-01"
[659] "2058-11-01" "2058-12-01" "2059-01-01" "2059-02-01" "2059-03-01" "2059-04-01" "2059-05-01"
[666] "2059-06-01" "2059-07-01" "2059-08-01" "2059-09-01" "2059-10-01" "2059-11-01" "2059-12-01"
[673] "2060-01-01" "2060-02-01" "2060-03-01" "2060-04-01" "2060-05-01" "2060-06-01" "2060-07-01"
[680] "2060-08-01" "2060-09-01" "2060-10-01" "2060-11-01" "2060-12-01" "2061-01-01" "2061-02-01"
[687] "2061-03-01" "2061-04-01" "2061-05-01" "2061-06-01" "2061-07-01" "2061-08-01" "2061-09-01"
[694] "2061-10-01" "2061-11-01" "2061-12-01" "2062-01-01" "2062-02-01" "2062-03-01" "2062-04-01"
[701] "2062-05-01" "2062-06-01" "2062-07-01" "2062-08-01" "2062-09-01" "2062-10-01" "2062-11-01"
[708] "2062-12-01" "2063-01-01" "2063-02-01" "2063-03-01" "2063-04-01" "2063-05-01" "2063-06-01"
[715] "2063-07-01" "2063-08-01" "2063-09-01" "2063-10-01" "2063-11-01" "2063-12-01" "2064-01-01"
[722] "2064-02-01" "2064-03-01" "2064-04-01" "2064-05-01" "2064-06-01" "2064-07-01" "2064-08-01"
[729] "2064-09-01" "2064-10-01" "2064-11-01" "2064-12-01" "2065-01-01" "2065-02-01" "2065-03-01"
[736] "2065-04-01" "2065-05-01" "2065-06-01" "2065-07-01" "2065-08-01" "2065-09-01" "2065-10-01"
[743] "2065-11-01" "2065-12-01" "2066-01-01" "2066-02-01" "2066-03-01" "2066-04-01" "2066-05-01"
[750] "2066-06-01" "2066-07-01" "2066-08-01" "2066-09-01" "2066-10-01" "2066-11-01" "2066-12-01"
[757] "2067-01-01" "2067-02-01" "2067-03-01" "2067-04-01" "2067-05-01" "2067-06-01" "2067-07-01"
[764] "2067-08-01" "2067-09-01" "2067-10-01" "2067-11-01" "2067-12-01" "2068-01-01" "2068-02-01"
[771] "2068-03-01" "2068-04-01" "2068-05-01" "2068-06-01" "2068-07-01" "2068-08-01" "2068-09-01"
[778] "2068-10-01" "2068-11-01" "2068-12-01" "2069-01-01" "2069-02-01" "2069-03-01" "2069-04-01"
[785] "2069-05-01" "2069-06-01" "2069-07-01" "2069-08-01" "2069-09-01" "2069-10-01" "2069-11-01"
[792] "2069-12-01" "2070-01-01" "2070-02-01" "2070-03-01" "2070-04-01" "2070-05-01" "2070-06-01"
[799] "2070-07-01" "2070-08-01" "2070-09-01" "2070-10-01" "2070-11-01" "2070-12-01" "2071-01-01"
[806] "2071-02-01" "2071-03-01" "2071-04-01" "2071-05-01" "2071-06-01" "2071-07-01" "2071-08-01"
[813] "2071-09-01" "2071-10-01" "2071-11-01" "2071-12-01" "2072-01-01" "2072-02-01" "2072-03-01"
[820] "2072-04-01" "2072-05-01" "2072-06-01" "2072-07-01" "2072-08-01" "2072-09-01" "2072-10-01"
[827] "2072-11-01" "2072-12-01" "2073-01-01" "2073-02-01" "2073-03-01" "2073-04-01" "2073-05-01"
[834] "2073-06-01" "2073-07-01" "2073-08-01" "2073-09-01" "2073-10-01" "2073-11-01" "2073-12-01"
[841] "2074-01-01" "2074-02-01" "2074-03-01" "2074-04-01" "2074-05-01" "2074-06-01" "2074-07-01"
[848] "2074-08-01" "2074-09-01" "2074-10-01" "2074-11-01" "2074-12-01" "2075-01-01" "2075-02-01"
[855] "2075-03-01" "2075-04-01" "2075-05-01" "2075-06-01" "2075-07-01" "2075-08-01" "2075-09-01"
[862] "2075-10-01" "2075-11-01" "2075-12-01" "2076-01-01" "2076-02-01" "2076-03-01" "2076-04-01"
[869] "2076-05-01" "2076-06-01" "2076-07-01" "2076-08-01" "2076-09-01" "2076-10-01" "2076-11-01"
[876] "2076-12-01" "2077-01-01" "2077-02-01" "2077-03-01" "2077-04-01" "2077-05-01" "2077-06-01"
[883] "2077-07-01" "2077-08-01" "2077-09-01" "2077-10-01" "2077-11-01" "2077-12-01" "2078-01-01"
[890] "2078-02-01" "2078-03-01" "2078-04-01" "2078-05-01" "2078-06-01" "2078-07-01" "2078-08-01"
[897] "2078-09-01" "2078-10-01" "2078-11-01" "2078-12-01" "2079-01-01" "2079-02-01" "2079-03-01"
[904] "2079-04-01" "2079-05-01" "2079-06-01" "2079-07-01" "2079-08-01" "2079-09-01" "2079-10-01"
[911] "2079-11-01" "2079-12-01" "2080-01-01" "2080-02-01" "2080-03-01" "2080-04-01" "2080-05-01"
[918] "2080-06-01" "2080-07-01" "2080-08-01" "2080-09-01" "2080-10-01" "2080-11-01" "2080-12-01"
[925] "2081-01-01" "2081-02-01" "2081-03-01" "2081-04-01" "2081-05-01" "2081-06-01" "2081-07-01"
[932] "2081-08-01" "2081-09-01" "2081-10-01" "2081-11-01" "2081-12-01" "2082-01-01" "2082-02-01"
[939] "2082-03-01" "2082-04-01" "2082-05-01" "2082-06-01" "2082-07-01" "2082-08-01" "2082-09-01"
[946] "2082-10-01" "2082-11-01" "2082-12-01" "2083-01-01" "2083-02-01" "2083-03-01" "2083-04-01"
[953] "2083-05-01" "2083-06-01" "2083-07-01" "2083-08-01" "2083-09-01" "2083-10-01" "2083-11-01"
[960] "2083-12-01" "2084-01-01" "2084-02-01" "2084-03-01" "2084-04-01" "2084-05-01" "2084-06-01"
[967] "2084-07-01" "2084-08-01" "2084-09-01" "2084-10-01" "2084-11-01" "2084-12-01" "2085-01-01"
[974] "2085-02-01" "2085-03-01" "2085-04-01" "2085-05-01" "2085-06-01" "2085-07-01" "2085-08-01"
[981] "2085-09-01" "2085-10-01" "2085-11-01" "2085-12-01" "2086-01-01" "2086-02-01" "2086-03-01"
[988] "2086-04-01" "2086-05-01" "2086-06-01" "2086-07-01" "2086-08-01" "2086-09-01" "2086-10-01"
[995] "2086-11-01" "2086-12-01" "2087-01-01" "2087-02-01" "2087-03-01" "2087-04-01"
[ reached getOption("max.print") -- omitted 164 entries ]
time <-as.Date(time)
out <- as.matrix(YF45_out_response) ## convert from coda to matrix
code from AMW for splitting up this matrix
x_cols <- grep('^x', colnames(out))
obs_cols <- grep('^OBS', colnames(out))
out_x <- out[,x_cols]
out_obs <- out[,obs_cols]
rm(out)
out_x <- t(out_x)
out_x <- as.data.frame(out_x)
out_x2 <- out_x |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*x\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
timesteps <- as.data.frame(time)
timesteps$ind <- as.character(seq(from = 1, to = nrow(timesteps), by = 1))
colnames(timesteps) <- c('time', 'timestep')
ponds <- as.data.frame(colnames(YF_WaterTemp45[-1]))
ponds$ind <- as.character(seq(from = 1, to = nrow(ponds), by = 1))
colnames(ponds) <- c('pond_name', 'pond')
out_x_mapped <- out_x2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pX <- out_x_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
out_obs <- t(out_obs)
out_obs <- as.data.frame(out_obs)
out_obs2 <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
out_obs_mapped <- out_obs2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pOBS <- out_obs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
cowplot::plot_grid(pX, pOBS)
threshold_probs <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(n_iter = dplyr::n(),
n_thresh = sum(val > 20),
prob = (n_thresh / n_iter) * 100)
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
threshold_probs_mapped <- threshold_probs |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond')
threshold_probs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = prob)) +
ggplot2::facet_wrap(~pond_name)
Save these ouputs using a .RDataFile type
save.image(file = "Corr_PooledMod_YF45.RData")
#Currently this has all the information in it. Eventually could limit this to be just the relavent information for the remaining code below
Load in this file here as long as no code changes are needed above
load(file = "Corr_PooledMod_YF45.RData")
Creating some nice plots of these
# set up the colors here for the plots below so that each pond has the same color each time
Mod_YF45 <- out_obs_mapped
str(Mod_YF45)
gropd_df [12,804 × 6] (S3: grouped_df/tbl_df/tbl/data.frame)
$ pond : chr [1:12804] "1" "1" "1" "1" ...
$ low : Named num [1:12804] -6.02 6.62 -1.38 1.07 4.39 ...
..- attr(*, "names")= chr [1:12804] "2.5%" "2.5%" "2.5%" "2.5%" ...
$ med : num [1:12804] 0.173 6.62 2.398 4.825 7.657 ...
$ high : Named num [1:12804] 6.4 6.62 5.98 8.39 10.77 ...
..- attr(*, "names")= chr [1:12804] "97.5%" "97.5%" "97.5%" "97.5%" ...
$ time : Date[1:12804], format: "2012-01-01" "2012-10-01" "2012-04-01" ...
$ pond_name: chr [1:12804] "MP1" "MP1" "MP1" "MP1" ...
- attr(*, "groups")= tibble [11 × 2] (S3: tbl_df/tbl/data.frame)
..$ pond : chr [1:11] "1" "10" "11" "2" ...
..$ .rows: list<int> [1:11]
.. ..$ : int [1:1164] 1 2 3 4 5 6 7 8 9 10 ...
.. ..$ : int [1:1164] 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 ...
.. ..$ : int [1:1164] 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 ...
.. ..$ : int [1:1164] 3493 3494 3495 3496 3497 3498 3499 3500 3501 3502 ...
.. ..$ : int [1:1164] 4657 4658 4659 4660 4661 4662 4663 4664 4665 4666 ...
.. ..$ : int [1:1164] 5821 5822 5823 5824 5825 5826 5827 5828 5829 5830 ...
.. ..$ : int [1:1164] 6985 6986 6987 6988 6989 6990 6991 6992 6993 6994 ...
.. ..$ : int [1:1164] 8149 8150 8151 8152 8153 8154 8155 8156 8157 8158 ...
.. ..$ : int [1:1164] 9313 9314 9315 9316 9317 9318 9319 9320 9321 9322 ...
.. ..$ : int [1:1164] 10477 10478 10479 10480 10481 10482 10483 10484 10485 10486 ...
.. ..$ : int [1:1164] 11641 11642 11643 11644 11645 11646 11647 11648 11649 11650 ...
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
write.csv(Mod_YF45, "Mod_YF45.csv")
# Define the custom color palette
YF_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "lightblue", "lightgreen", "#6A3D9A")
# Ensure that YF_colors has the same length as the number of unique ponds
unique_ponds <- unique(Mod_YF45$pond_name)
if (length(YF_colors) != length(unique_ponds)) {
stop("The number of colors in YF_colors does not match the number of unique ponds.")
}
# Create a named vector for YF_colors with pond names
color_mapping <- setNames(YF_colors, unique_ponds)
# Plotting
Ponds_YF45 <- ggplot(data = Mod_YF45, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Add the ribbon with colors
geom_line(aes(y = med, color = pond_name)) + # Ensure lines are visible and colored by pond
geom_smooth(aes(y = med, color = "grey"), size = 1) + # Ensure lines are visible and colored by pond
geom_vline(xintercept = as.Date("2020-12-01"), linetype = "dashed", color = "red") + # Add vertical dashed red line
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
xlab("Year") +
ylab("Water Temperature (°C)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Ponds_YF45
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndPonds_YF45.png", plot = Ponds_YF45, width = 8, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
# using the same aestetics as above
Trends_YF45 <- ggplot(data = Mod_YF45, aes(x = time)) +
geom_smooth(aes(y = med, color = pond_name, fill = pond_name), alpha = 0.15, size = 1) + # Smooth lines with colors by pond_name
scale_color_manual(values = color_mapping) +
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
legend.position = "right",
legend.text = element_text(size = 16), # Change the legend text size here
legend.title = element_text(size = 20))
Trends_YF45
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndTrends_YF45.png", plot = Trends_YF45, width = 9, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
probs_YF45 <- threshold_probs_mapped
write.csv(probs_YF45, "probs_YF45.csv")
Probs20_YF45 <- ggplot(data = probs_YF45, aes(x = time)) +
geom_line(aes(y = prob, color = pond_name)) + # Ensure lines are visible and colored by pond
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
xlab("Year") +
ylab("Probability") +
xlim(as.Date(c("2020-12-01", "2099-12-01"))) +
ylim(0,50) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Probs20_YF45
Warning: Removed 2365 rows containing missing values or values outside the scale range (`geom_line()`).
ggsave("Corr_Probs20_YF45.png", plot = Probs20_YF45, width = 10, height = 6)
Warning: Removed 2365 rows containing missing values or values outside the scale range (`geom_line()`).
Mod_YF45
ModObs <- Mod_YF45 %>%
filter(time < "2021-01-01")
range(ModObs$time)
[1] "2012-01-01" "2020-12-01"
ErrorPlot_YF45 <- ggplot(data = ModObs, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Adjust fill colors
geom_line(aes(y = med), color = "black") + # Black lines for the median values
facet_wrap(~ pond_name) + # Facet by pond_name
scale_fill_manual(values = color_mapping) + # Apply custom color palette for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(
plot.title = element_text(hjust = 0.5),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none" # Angle x-axis labels if needed
)
ErrorPlot_YF45
ggsave("Corr_Error_YF45.png", plot = ErrorPlot_YF45, width = 9, height = 6)
Clear out the environment and reload the basics before starting the next section
rm(list = ls())
library(tidyverse)
library(rjags)
library(ggplot2)
set.seed(1)
Read in the data here, which is an output of XX code file
Then separate the data into air and water temperature files
# Air and water temperature data: YF loccation, 4.5 scenario
YF_Temps <- read.csv("YF_85SNAPForecast.csv", header=TRUE)
# Air temperature variables
YF_AirTemp85 <- YF_Temps %>%
select(date, Airport, Air_MonthAvg, std, precision, Region)
# Water temperature variables
YF_WaterTemp85 <- YF_Temps %>%
select(date, MP1, MP3, MP5, MP8, PL1, PL2, PL3, UBP1, UBP2, UBP3, UBP4)
# Set up a string of all the dates
time <- YF_WaterTemp85$date
length(time)
[1] 1164
Setting up the random walk model
dlm_pooled <- "
model{
#### Priors
for(p in 1:np){
x[1,p] ~ dnorm(x_ic, tau_ic) # Initial condition of water temperature
}
tau_obs ~ dgamma(a_obs, r_obs) # Prior on observation error
tau_add ~ dgamma(a_add, r_add) # Prior on process error
#### Fixed Effects
beta ~ dmnorm(mu_beta, tau_beta) # Prior on beta coefficients
for(p in 1:np){
int[p] ~ dnorm(mu_int, tau_int) # Prior on pond-specific intercepts
}
#### Data Model
for(t in 1:n){ # loop over all time steps
for(p in 1:np){
OBS[t,p] ~ dnorm(x[t,p], tau_obs) # Observed water temperature is drawn from latent air temperature with observation uncertainty
}
Xf[t] ~ dnorm(muAirTemp[t], tauAirTemp[t]) # Latent air temperature is drawn from mean and precision of forecasated air temperature
}
#### Process Model
for(t in 2:n){ # loop over all time steps except teh first (we defined ic above)
for(p in 1:np){
mu[t,p] <- x[t-1,p] + int[p] + beta[1] * x[t-1,p] + beta[2] * Xf[t] # Mean water temperature is a function of the previous time step and current air temperature
x[t,p] ~ dnorm(mu[t,p], tau_add) # Latent water temperature is drawn from mean water temperature with process uncertainty
}
}
}
"
Define the data and priors for the model
# Empty list
data <- list()
# Water temperature observations
data$OBS <- dplyr::select(YF_WaterTemp85, -date)
# Number of time steps
data$n <- nrow(data$OBS)
# Number of ponds
data$np <- ncol(data$OBS)
# Initial water temperature mean
data$x_ic <- 0.1
# Initial water temperature precision
data$tau_ic = 0.1
# Prior parameters for observation and process uncertainty
data$a_obs = 1
data$r_obs = 1
data$a_add = 1
data$r_add = 1
# Prior parameters for beta coefficients
data$mu_beta <- c(0, 0)
data$tau_beta <- diag(x = c(0.001, 0.001), nrow = 2, ncol = 2)
# Prior parameters for intercept
data$mu_int <- 0
data$tau_int <- 0.001
# Mean air temperature estimate
data$muAirTemp <- YF_AirTemp85$Air_MonthAvg
# Air temperature precision
data$tauAirTemp <- YF_AirTemp85$precision
Create JAGS model with 3 chains
jm <- jags.model(file = textConnection(dlm_pooled), data = data, n.chains = 3)
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 456
Unobserved stochastic nodes: 26330
Total graph size: 55883
Initializing model
Posterior samples of parameters
YF85_out_params <- coda.samples(model = jm,
variable.names = c('beta', 'int',
'tau_add', 'tau_obs'),
n.iter = 50000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
Posterior samples of response variables
YF85_out_response <- coda.samples(model = jm,
variable.names = c('x', 'OBS'),
n.iter = 100000, thin = 25)
|
| | 0%
|
|* | 2%
|
|** | 4%
|
|*** | 6%
|
|**** | 8%
|
|***** | 10%
|
|****** | 12%
|
|******* | 14%
|
|******** | 16%
|
|********* | 18%
|
|********** | 20%
|
|*********** | 22%
|
|************ | 24%
|
|************* | 26%
|
|************** | 28%
|
|*************** | 30%
|
|**************** | 32%
|
|***************** | 34%
|
|****************** | 36%
|
|******************* | 38%
|
|******************** | 40%
|
|********************* | 42%
|
|********************** | 44%
|
|*********************** | 46%
|
|************************ | 48%
|
|************************* | 50%
|
|************************** | 52%
|
|*************************** | 54%
|
|**************************** | 56%
|
|***************************** | 58%
|
|****************************** | 60%
|
|******************************* | 62%
|
|******************************** | 64%
|
|********************************* | 66%
|
|********************************** | 68%
|
|*********************************** | 70%
|
|************************************ | 72%
|
|************************************* | 74%
|
|************************************** | 76%
|
|*************************************** | 78%
|
|**************************************** | 80%
|
|***************************************** | 82%
|
|****************************************** | 84%
|
|******************************************* | 86%
|
|******************************************** | 88%
|
|********************************************* | 90%
|
|********************************************** | 92%
|
|*********************************************** | 94%
|
|************************************************ | 96%
|
|************************************************* | 98%
|
|**************************************************| 100%
A couple quick checks here
# plot(YF85_out_params)
gelman.diag(YF85_out_params, confidence = 0.99)
Potential scale reduction factors:
Point est. Upper C.I.
beta[1] 1.00 1.01
beta[2] 1.00 1.01
int[1] 1.01 1.03
int[2] 1.01 1.05
int[3] 1.01 1.04
int[4] 1.01 1.04
int[5] 1.01 1.04
int[6] 1.01 1.05
int[7] 1.01 1.05
int[8] 1.01 1.04
int[9] 1.01 1.05
int[10] 1.01 1.05
int[11] 1.01 1.04
tau_add 1.00 1.00
tau_obs 1.00 1.00
Multivariate psrf
1.01
Visualize the output by just looking at the 95% Credible interval of the time-series of X’s and compare that to the observed Y’s
Transform the samples back from the log domain to the linear domain
time ## adjust to zoom in and out
[1] "2012-01-01" "2012-02-01" "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01"
[8] "2012-08-01" "2012-09-01" "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01"
[15] "2013-03-01" "2013-04-01" "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01"
[22] "2013-10-01" "2013-11-01" "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01"
[29] "2014-05-01" "2014-06-01" "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01"
[36] "2014-12-01" "2015-01-01" "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01"
[43] "2015-07-01" "2015-08-01" "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01"
[50] "2016-02-01" "2016-03-01" "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01"
[57] "2016-09-01" "2016-10-01" "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01"
[64] "2017-04-01" "2017-05-01" "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01"
[71] "2017-11-01" "2017-12-01" "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01"
[78] "2018-06-01" "2018-07-01" "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01"
[85] "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01"
[92] "2019-08-01" "2019-09-01" "2019-10-01" "2019-11-01" "2019-12-01" "2012-01-01" "2012-02-01"
[99] "2012-03-01" "2012-04-01" "2012-05-01" "2012-06-01" "2012-07-01" "2012-08-01" "2012-09-01"
[106] "2012-10-01" "2012-11-01" "2012-12-01" "2013-01-01" "2013-02-01" "2013-03-01" "2013-04-01"
[113] "2013-05-01" "2013-06-01" "2013-07-01" "2013-08-01" "2013-09-01" "2013-10-01" "2013-11-01"
[120] "2013-12-01" "2014-01-01" "2014-02-01" "2014-03-01" "2014-04-01" "2014-05-01" "2014-06-01"
[127] "2014-07-01" "2014-08-01" "2014-09-01" "2014-10-01" "2014-11-01" "2014-12-01" "2015-01-01"
[134] "2015-02-01" "2015-03-01" "2015-04-01" "2015-05-01" "2015-06-01" "2015-07-01" "2015-08-01"
[141] "2015-09-01" "2015-10-01" "2015-11-01" "2015-12-01" "2016-01-01" "2016-02-01" "2016-03-01"
[148] "2016-04-01" "2016-05-01" "2016-06-01" "2016-07-01" "2016-08-01" "2016-09-01" "2016-10-01"
[155] "2016-11-01" "2016-12-01" "2017-01-01" "2017-02-01" "2017-03-01" "2017-04-01" "2017-05-01"
[162] "2017-06-01" "2017-07-01" "2017-08-01" "2017-09-01" "2017-10-01" "2017-11-01" "2017-12-01"
[169] "2018-01-01" "2018-02-01" "2018-03-01" "2018-04-01" "2018-05-01" "2018-06-01" "2018-07-01"
[176] "2018-08-01" "2018-09-01" "2018-10-01" "2018-11-01" "2018-12-01" "2019-01-01" "2019-02-01"
[183] "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01" "2019-07-01" "2019-08-01" "2019-09-01"
[190] "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01" "2020-03-01" "2020-04-01"
[197] "2020-05-01" "2020-06-01" "2020-07-01" "2020-08-01" "2020-09-01" "2020-10-01" "2020-11-01"
[204] "2020-12-01" "2021-01-01" "2021-02-01" "2021-03-01" "2021-04-01" "2021-05-01" "2021-06-01"
[211] "2021-07-01" "2021-08-01" "2021-09-01" "2021-10-01" "2021-11-01" "2021-12-01" "2022-01-01"
[218] "2022-02-01" "2022-03-01" "2022-04-01" "2022-05-01" "2022-06-01" "2022-07-01" "2022-08-01"
[225] "2022-09-01" "2022-10-01" "2022-11-01" "2022-12-01" "2023-01-01" "2023-02-01" "2023-03-01"
[232] "2023-04-01" "2023-05-01" "2023-06-01" "2023-07-01" "2023-08-01" "2023-09-01" "2023-10-01"
[239] "2023-11-01" "2023-12-01" "2024-01-01" "2024-02-01" "2024-03-01" "2024-04-01" "2024-05-01"
[246] "2024-06-01" "2024-07-01" "2024-08-01" "2024-09-01" "2024-10-01" "2024-11-01" "2024-12-01"
[253] "2025-01-01" "2025-02-01" "2025-03-01" "2025-04-01" "2025-05-01" "2025-06-01" "2025-07-01"
[260] "2025-08-01" "2025-09-01" "2025-10-01" "2025-11-01" "2025-12-01" "2026-01-01" "2026-02-01"
[267] "2026-03-01" "2026-04-01" "2026-05-01" "2026-06-01" "2026-07-01" "2026-08-01" "2026-09-01"
[274] "2026-10-01" "2026-11-01" "2026-12-01" "2027-01-01" "2027-02-01" "2027-03-01" "2027-04-01"
[281] "2027-05-01" "2027-06-01" "2027-07-01" "2027-08-01" "2027-09-01" "2027-10-01" "2027-11-01"
[288] "2027-12-01" "2028-01-01" "2028-02-01" "2028-03-01" "2028-04-01" "2028-05-01" "2028-06-01"
[295] "2028-07-01" "2028-08-01" "2028-09-01" "2028-10-01" "2028-11-01" "2028-12-01" "2029-01-01"
[302] "2029-02-01" "2029-03-01" "2029-04-01" "2029-05-01" "2029-06-01" "2029-07-01" "2029-08-01"
[309] "2029-09-01" "2029-10-01" "2029-11-01" "2029-12-01" "2030-01-01" "2030-02-01" "2030-03-01"
[316] "2030-04-01" "2030-05-01" "2030-06-01" "2030-07-01" "2030-08-01" "2030-09-01" "2030-10-01"
[323] "2030-11-01" "2030-12-01" "2031-01-01" "2031-02-01" "2031-03-01" "2031-04-01" "2031-05-01"
[330] "2031-06-01" "2031-07-01" "2031-08-01" "2031-09-01" "2031-10-01" "2031-11-01" "2031-12-01"
[337] "2032-01-01" "2032-02-01" "2032-03-01" "2032-04-01" "2032-05-01" "2032-06-01" "2032-07-01"
[344] "2032-08-01" "2032-09-01" "2032-10-01" "2032-11-01" "2032-12-01" "2033-01-01" "2033-02-01"
[351] "2033-03-01" "2033-04-01" "2033-05-01" "2033-06-01" "2033-07-01" "2033-08-01" "2033-09-01"
[358] "2033-10-01" "2033-11-01" "2033-12-01" "2034-01-01" "2034-02-01" "2034-03-01" "2034-04-01"
[365] "2034-05-01" "2034-06-01" "2034-07-01" "2034-08-01" "2034-09-01" "2034-10-01" "2034-11-01"
[372] "2034-12-01" "2035-01-01" "2035-02-01" "2035-03-01" "2035-04-01" "2035-05-01" "2035-06-01"
[379] "2035-07-01" "2035-08-01" "2035-09-01" "2035-10-01" "2035-11-01" "2035-12-01" "2036-01-01"
[386] "2036-02-01" "2036-03-01" "2036-04-01" "2036-05-01" "2036-06-01" "2036-07-01" "2036-08-01"
[393] "2036-09-01" "2036-10-01" "2036-11-01" "2036-12-01" "2037-01-01" "2037-02-01" "2037-03-01"
[400] "2037-04-01" "2037-05-01" "2037-06-01" "2037-07-01" "2037-08-01" "2037-09-01" "2037-10-01"
[407] "2037-11-01" "2037-12-01" "2038-01-01" "2038-02-01" "2038-03-01" "2038-04-01" "2038-05-01"
[414] "2038-06-01" "2038-07-01" "2038-08-01" "2038-09-01" "2038-10-01" "2038-11-01" "2038-12-01"
[421] "2039-01-01" "2039-02-01" "2039-03-01" "2039-04-01" "2039-05-01" "2039-06-01" "2039-07-01"
[428] "2039-08-01" "2039-09-01" "2039-10-01" "2039-11-01" "2039-12-01" "2040-01-01" "2040-02-01"
[435] "2040-03-01" "2040-04-01" "2040-05-01" "2040-06-01" "2040-07-01" "2040-08-01" "2040-09-01"
[442] "2040-10-01" "2040-11-01" "2040-12-01" "2041-01-01" "2041-02-01" "2041-03-01" "2041-04-01"
[449] "2041-05-01" "2041-06-01" "2041-07-01" "2041-08-01" "2041-09-01" "2041-10-01" "2041-11-01"
[456] "2041-12-01" "2042-01-01" "2042-02-01" "2042-03-01" "2042-04-01" "2042-05-01" "2042-06-01"
[463] "2042-07-01" "2042-08-01" "2042-09-01" "2042-10-01" "2042-11-01" "2042-12-01" "2043-01-01"
[470] "2043-02-01" "2043-03-01" "2043-04-01" "2043-05-01" "2043-06-01" "2043-07-01" "2043-08-01"
[477] "2043-09-01" "2043-10-01" "2043-11-01" "2043-12-01" "2044-01-01" "2044-02-01" "2044-03-01"
[484] "2044-04-01" "2044-05-01" "2044-06-01" "2044-07-01" "2044-08-01" "2044-09-01" "2044-10-01"
[491] "2044-11-01" "2044-12-01" "2045-01-01" "2045-02-01" "2045-03-01" "2045-04-01" "2045-05-01"
[498] "2045-06-01" "2045-07-01" "2045-08-01" "2045-09-01" "2045-10-01" "2045-11-01" "2045-12-01"
[505] "2046-01-01" "2046-02-01" "2046-03-01" "2046-04-01" "2046-05-01" "2046-06-01" "2046-07-01"
[512] "2046-08-01" "2046-09-01" "2046-10-01" "2046-11-01" "2046-12-01" "2047-01-01" "2047-02-01"
[519] "2047-03-01" "2047-04-01" "2047-05-01" "2047-06-01" "2047-07-01" "2047-08-01" "2047-09-01"
[526] "2047-10-01" "2047-11-01" "2047-12-01" "2048-01-01" "2048-02-01" "2048-03-01" "2048-04-01"
[533] "2048-05-01" "2048-06-01" "2048-07-01" "2048-08-01" "2048-09-01" "2048-10-01" "2048-11-01"
[540] "2048-12-01" "2049-01-01" "2049-02-01" "2049-03-01" "2049-04-01" "2049-05-01" "2049-06-01"
[547] "2049-07-01" "2049-08-01" "2049-09-01" "2049-10-01" "2049-11-01" "2049-12-01" "2050-01-01"
[554] "2050-02-01" "2050-03-01" "2050-04-01" "2050-05-01" "2050-06-01" "2050-07-01" "2050-08-01"
[561] "2050-09-01" "2050-10-01" "2050-11-01" "2050-12-01" "2051-01-01" "2051-02-01" "2051-03-01"
[568] "2051-04-01" "2051-05-01" "2051-06-01" "2051-07-01" "2051-08-01" "2051-09-01" "2051-10-01"
[575] "2051-11-01" "2051-12-01" "2052-01-01" "2052-02-01" "2052-03-01" "2052-04-01" "2052-05-01"
[582] "2052-06-01" "2052-07-01" "2052-08-01" "2052-09-01" "2052-10-01" "2052-11-01" "2052-12-01"
[589] "2053-01-01" "2053-02-01" "2053-03-01" "2053-04-01" "2053-05-01" "2053-06-01" "2053-07-01"
[596] "2053-08-01" "2053-09-01" "2053-10-01" "2053-11-01" "2053-12-01" "2054-01-01" "2054-02-01"
[603] "2054-03-01" "2054-04-01" "2054-05-01" "2054-06-01" "2054-07-01" "2054-08-01" "2054-09-01"
[610] "2054-10-01" "2054-11-01" "2054-12-01" "2055-01-01" "2055-02-01" "2055-03-01" "2055-04-01"
[617] "2055-05-01" "2055-06-01" "2055-07-01" "2055-08-01" "2055-09-01" "2055-10-01" "2055-11-01"
[624] "2055-12-01" "2056-01-01" "2056-02-01" "2056-03-01" "2056-04-01" "2056-05-01" "2056-06-01"
[631] "2056-07-01" "2056-08-01" "2056-09-01" "2056-10-01" "2056-11-01" "2056-12-01" "2057-01-01"
[638] "2057-02-01" "2057-03-01" "2057-04-01" "2057-05-01" "2057-06-01" "2057-07-01" "2057-08-01"
[645] "2057-09-01" "2057-10-01" "2057-11-01" "2057-12-01" "2058-01-01" "2058-02-01" "2058-03-01"
[652] "2058-04-01" "2058-05-01" "2058-06-01" "2058-07-01" "2058-08-01" "2058-09-01" "2058-10-01"
[659] "2058-11-01" "2058-12-01" "2059-01-01" "2059-02-01" "2059-03-01" "2059-04-01" "2059-05-01"
[666] "2059-06-01" "2059-07-01" "2059-08-01" "2059-09-01" "2059-10-01" "2059-11-01" "2059-12-01"
[673] "2060-01-01" "2060-02-01" "2060-03-01" "2060-04-01" "2060-05-01" "2060-06-01" "2060-07-01"
[680] "2060-08-01" "2060-09-01" "2060-10-01" "2060-11-01" "2060-12-01" "2061-01-01" "2061-02-01"
[687] "2061-03-01" "2061-04-01" "2061-05-01" "2061-06-01" "2061-07-01" "2061-08-01" "2061-09-01"
[694] "2061-10-01" "2061-11-01" "2061-12-01" "2062-01-01" "2062-02-01" "2062-03-01" "2062-04-01"
[701] "2062-05-01" "2062-06-01" "2062-07-01" "2062-08-01" "2062-09-01" "2062-10-01" "2062-11-01"
[708] "2062-12-01" "2063-01-01" "2063-02-01" "2063-03-01" "2063-04-01" "2063-05-01" "2063-06-01"
[715] "2063-07-01" "2063-08-01" "2063-09-01" "2063-10-01" "2063-11-01" "2063-12-01" "2064-01-01"
[722] "2064-02-01" "2064-03-01" "2064-04-01" "2064-05-01" "2064-06-01" "2064-07-01" "2064-08-01"
[729] "2064-09-01" "2064-10-01" "2064-11-01" "2064-12-01" "2065-01-01" "2065-02-01" "2065-03-01"
[736] "2065-04-01" "2065-05-01" "2065-06-01" "2065-07-01" "2065-08-01" "2065-09-01" "2065-10-01"
[743] "2065-11-01" "2065-12-01" "2066-01-01" "2066-02-01" "2066-03-01" "2066-04-01" "2066-05-01"
[750] "2066-06-01" "2066-07-01" "2066-08-01" "2066-09-01" "2066-10-01" "2066-11-01" "2066-12-01"
[757] "2067-01-01" "2067-02-01" "2067-03-01" "2067-04-01" "2067-05-01" "2067-06-01" "2067-07-01"
[764] "2067-08-01" "2067-09-01" "2067-10-01" "2067-11-01" "2067-12-01" "2068-01-01" "2068-02-01"
[771] "2068-03-01" "2068-04-01" "2068-05-01" "2068-06-01" "2068-07-01" "2068-08-01" "2068-09-01"
[778] "2068-10-01" "2068-11-01" "2068-12-01" "2069-01-01" "2069-02-01" "2069-03-01" "2069-04-01"
[785] "2069-05-01" "2069-06-01" "2069-07-01" "2069-08-01" "2069-09-01" "2069-10-01" "2069-11-01"
[792] "2069-12-01" "2070-01-01" "2070-02-01" "2070-03-01" "2070-04-01" "2070-05-01" "2070-06-01"
[799] "2070-07-01" "2070-08-01" "2070-09-01" "2070-10-01" "2070-11-01" "2070-12-01" "2071-01-01"
[806] "2071-02-01" "2071-03-01" "2071-04-01" "2071-05-01" "2071-06-01" "2071-07-01" "2071-08-01"
[813] "2071-09-01" "2071-10-01" "2071-11-01" "2071-12-01" "2072-01-01" "2072-02-01" "2072-03-01"
[820] "2072-04-01" "2072-05-01" "2072-06-01" "2072-07-01" "2072-08-01" "2072-09-01" "2072-10-01"
[827] "2072-11-01" "2072-12-01" "2073-01-01" "2073-02-01" "2073-03-01" "2073-04-01" "2073-05-01"
[834] "2073-06-01" "2073-07-01" "2073-08-01" "2073-09-01" "2073-10-01" "2073-11-01" "2073-12-01"
[841] "2074-01-01" "2074-02-01" "2074-03-01" "2074-04-01" "2074-05-01" "2074-06-01" "2074-07-01"
[848] "2074-08-01" "2074-09-01" "2074-10-01" "2074-11-01" "2074-12-01" "2075-01-01" "2075-02-01"
[855] "2075-03-01" "2075-04-01" "2075-05-01" "2075-06-01" "2075-07-01" "2075-08-01" "2075-09-01"
[862] "2075-10-01" "2075-11-01" "2075-12-01" "2076-01-01" "2076-02-01" "2076-03-01" "2076-04-01"
[869] "2076-05-01" "2076-06-01" "2076-07-01" "2076-08-01" "2076-09-01" "2076-10-01" "2076-11-01"
[876] "2076-12-01" "2077-01-01" "2077-02-01" "2077-03-01" "2077-04-01" "2077-05-01" "2077-06-01"
[883] "2077-07-01" "2077-08-01" "2077-09-01" "2077-10-01" "2077-11-01" "2077-12-01" "2078-01-01"
[890] "2078-02-01" "2078-03-01" "2078-04-01" "2078-05-01" "2078-06-01" "2078-07-01" "2078-08-01"
[897] "2078-09-01" "2078-10-01" "2078-11-01" "2078-12-01" "2079-01-01" "2079-02-01" "2079-03-01"
[904] "2079-04-01" "2079-05-01" "2079-06-01" "2079-07-01" "2079-08-01" "2079-09-01" "2079-10-01"
[911] "2079-11-01" "2079-12-01" "2080-01-01" "2080-02-01" "2080-03-01" "2080-04-01" "2080-05-01"
[918] "2080-06-01" "2080-07-01" "2080-08-01" "2080-09-01" "2080-10-01" "2080-11-01" "2080-12-01"
[925] "2081-01-01" "2081-02-01" "2081-03-01" "2081-04-01" "2081-05-01" "2081-06-01" "2081-07-01"
[932] "2081-08-01" "2081-09-01" "2081-10-01" "2081-11-01" "2081-12-01" "2082-01-01" "2082-02-01"
[939] "2082-03-01" "2082-04-01" "2082-05-01" "2082-06-01" "2082-07-01" "2082-08-01" "2082-09-01"
[946] "2082-10-01" "2082-11-01" "2082-12-01" "2083-01-01" "2083-02-01" "2083-03-01" "2083-04-01"
[953] "2083-05-01" "2083-06-01" "2083-07-01" "2083-08-01" "2083-09-01" "2083-10-01" "2083-11-01"
[960] "2083-12-01" "2084-01-01" "2084-02-01" "2084-03-01" "2084-04-01" "2084-05-01" "2084-06-01"
[967] "2084-07-01" "2084-08-01" "2084-09-01" "2084-10-01" "2084-11-01" "2084-12-01" "2085-01-01"
[974] "2085-02-01" "2085-03-01" "2085-04-01" "2085-05-01" "2085-06-01" "2085-07-01" "2085-08-01"
[981] "2085-09-01" "2085-10-01" "2085-11-01" "2085-12-01" "2086-01-01" "2086-02-01" "2086-03-01"
[988] "2086-04-01" "2086-05-01" "2086-06-01" "2086-07-01" "2086-08-01" "2086-09-01" "2086-10-01"
[995] "2086-11-01" "2086-12-01" "2087-01-01" "2087-02-01" "2087-03-01" "2087-04-01"
[ reached getOption("max.print") -- omitted 164 entries ]
time <-as.Date(time)
out <- as.matrix(YF85_out_response) ## convert from coda to matrix
code from AMW for splitting up this matrix
x_cols <- grep('^x', colnames(out))
obs_cols <- grep('^OBS', colnames(out))
out_x <- out[,x_cols]
out_obs <- out[,obs_cols]
rm(out)
out_x <- t(out_x)
out_x <- as.data.frame(out_x)
out_x2 <- out_x |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*x\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
timesteps <- as.data.frame(time)
timesteps$ind <- as.character(seq(from = 1, to = nrow(timesteps), by = 1))
colnames(timesteps) <- c('time', 'timestep')
ponds <- as.data.frame(colnames(YF_WaterTemp85[-1]))
ponds$ind <- as.character(seq(from = 1, to = nrow(ponds), by = 1))
colnames(ponds) <- c('pond_name', 'pond')
out_x_mapped <- out_x2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pX <- out_x_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
out_obs <- t(out_obs)
out_obs <- as.data.frame(out_obs)
out_obs2 <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(low = quantile(val, probs = 0.025),
med = median(val),
high = quantile(val, probs = 0.975))
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
out_obs_mapped <- out_obs2 |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond') |>
dplyr::select(low, med, high, time, pond_name)
Adding missing grouping variables: `pond`
pOBS <- out_obs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = med)) +
ggplot2::geom_ribbon(ggplot2::aes(x = time, ymin = low, ymax = high)) +
ggplot2::facet_wrap(~pond_name)
cowplot::plot_grid(pX, pOBS)
threshold_probs <- out_obs |>
tibble::rownames_to_column(var = 'time_pond') |>
dplyr::mutate(timestep = sub(x = time_pond,
pattern = '.*OBS\\[',
replacement = ''),
timestep = sub(x = timestep,
pattern = ',.*',
replacement = ''),
pond = sub(x = time_pond,
pattern = '.*,',
replacement = ''),
pond = sub(x = pond,
pattern = '\\]',
replacement = '')) |>
tidyr::pivot_longer(-c(time_pond, timestep, pond),
names_to = 'iter', values_to = 'val') |>
dplyr::group_by(pond, timestep) |>
dplyr::summarize(n_iter = dplyr::n(),
n_thresh = sum(val > 20),
prob = (n_thresh / n_iter) * 100)
`summarise()` has grouped output by 'pond'. You can override using the `.groups` argument.
threshold_probs_mapped <- threshold_probs |>
dplyr::left_join(y = timesteps, by = 'timestep') |>
dplyr::left_join(y = ponds, by = 'pond')
threshold_probs_mapped |>
ggplot2::ggplot() +
ggplot2::geom_line(ggplot2::aes(x = time, y = prob)) +
ggplot2::facet_wrap(~pond_name)
Save these ouputs using a .RDataFile type
save.image(file = "Corr_PooledMod_YF85.RData")
rsession-arm64(14651) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
#Currently this has all the information in it. Eventually could limit this to be just the relavent information for the remaining code below
Load in this file here as long as no code changes are needed above
load(file = "Corr_PooledMod_YF85.RData")
Creating some nice plots of these
# set up the colors here for the plots below so that each pond has the same color each time
Mod_YF85 <- out_obs_mapped
str(Mod_YF85)
gropd_df [12,804 × 6] (S3: grouped_df/tbl_df/tbl/data.frame)
$ pond : chr [1:12804] "1" "1" "1" "1" ...
$ low : Named num [1:12804] -6.03 6.62 -2.05 2.66 5.44 ...
..- attr(*, "names")= chr [1:12804] "2.5%" "2.5%" "2.5%" "2.5%" ...
$ med : num [1:12804] 0.166 6.62 2.018 6.625 9.66 ...
$ high : Named num [1:12804] 6.43 6.62 5.94 10.46 13.96 ...
..- attr(*, "names")= chr [1:12804] "97.5%" "97.5%" "97.5%" "97.5%" ...
$ time : Date[1:12804], format: "2012-01-01" "2012-10-01" "2012-04-01" ...
$ pond_name: chr [1:12804] "MP1" "MP1" "MP1" "MP1" ...
- attr(*, "groups")= tibble [11 × 2] (S3: tbl_df/tbl/data.frame)
..$ pond : chr [1:11] "1" "10" "11" "2" ...
..$ .rows: list<int> [1:11]
.. ..$ : int [1:1164] 1 2 3 4 5 6 7 8 9 10 ...
.. ..$ : int [1:1164] 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 ...
.. ..$ : int [1:1164] 2329 2330 2331 2332 2333 2334 2335 2336 2337 2338 ...
.. ..$ : int [1:1164] 3493 3494 3495 3496 3497 3498 3499 3500 3501 3502 ...
.. ..$ : int [1:1164] 4657 4658 4659 4660 4661 4662 4663 4664 4665 4666 ...
.. ..$ : int [1:1164] 5821 5822 5823 5824 5825 5826 5827 5828 5829 5830 ...
.. ..$ : int [1:1164] 6985 6986 6987 6988 6989 6990 6991 6992 6993 6994 ...
.. ..$ : int [1:1164] 8149 8150 8151 8152 8153 8154 8155 8156 8157 8158 ...
.. ..$ : int [1:1164] 9313 9314 9315 9316 9317 9318 9319 9320 9321 9322 ...
.. ..$ : int [1:1164] 10477 10478 10479 10480 10481 10482 10483 10484 10485 10486 ...
.. ..$ : int [1:1164] 11641 11642 11643 11644 11645 11646 11647 11648 11649 11650 ...
.. ..@ ptype: int(0)
..- attr(*, ".drop")= logi TRUE
write.csv(Mod_YF85, "Mod_YF85.csv")
rsession-arm64(14768) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
# Define the custom color palette
YF_colors <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "lightblue", "lightgreen", "#6A3D9A")
# Ensure that YF_colors has the same length as the number of unique ponds
unique_ponds <- unique(Mod_YF85$pond_name)
if (length(YF_colors) != length(unique_ponds)) {
stop("The number of colors in YF_colors does not match the number of unique ponds.")
}
# Create a named vector for YF_colors with pond names
color_mapping <- setNames(YF_colors, unique_ponds)
# Plotting
Ponds_YF85 <- ggplot(data = Mod_YF85, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Add the ribbon with colors
geom_line(aes(y = med, color = pond_name)) + # Ensure lines are visible and colored by pond
geom_smooth(aes(y = med, color = "grey"), size = 1) + # Ensure lines are visible and colored by pond
geom_vline(xintercept = as.Date("2020-12-01"), linetype = "dashed", color = "red") + # Add vertical dashed red line
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
xlab("Year") +
ylab("Water Temperature (°C)") +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
rsession-arm64(14787) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
Ponds_YF85
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
ggsave("Corr_IndPonds_YF85.png", plot = Ponds_YF85, width = 8, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
Create plots of the trend by pond (geom smooth for each pond on the same plot for comparison)
# using the same aestetics as above
Trends_YF85 <- ggplot(data = Mod_YF85, aes(x = time)) +
geom_smooth(aes(y = med, color = pond_name, fill = pond_name), alpha = 0.15, size = 1) + # Smooth lines with colors by pond_name
scale_color_manual(values = color_mapping) +
scale_fill_manual(values = color_mapping) + # Use custom colors for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
legend.position = "right",
legend.text = element_text(size = 16), # Change the legend text size here
legend.title = element_text(size = 20))
Trends_YF85
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
rsession-arm64(14797) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
ggsave("Corr_IndTrends_YF85.png", plot = Trends_YF85, width = 9, height = 6)
`geom_smooth()` using method = 'gam' and formula = 'y ~ s(x, bs = "cs")'
probs_YF85 <- threshold_probs_mapped
write.csv(probs_YF85, "probs_YF85.csv")
Probs20_YF85 <- ggplot(data = probs_YF85, aes(x = time)) +
geom_line(aes(y = prob, color = pond_name)) + # Ensure lines are visible and colored by pond
facet_wrap(~ pond_name) + # Facet by pond_name
scale_color_manual(values = color_mapping) + # Use custom colors for lines
xlab("Year") +
ylab("Probability") +
xlim(as.Date(c("2020-12-01", "2099-12-01"))) +
ylim(0,50) +
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 24),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
Probs20_YF85
Warning: Removed 2365 rows containing missing values or values outside the scale range (`geom_line()`).
rsession-arm64(14803) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
ggsave("Corr_Probs20_YF85.png", plot = Probs20_YF85, width = 9, height = 6)
Warning: Removed 2365 rows containing missing values or values outside the scale range (`geom_line()`).
Mod_YF85
ModObs <- Mod_YF85 %>%
filter(time < "2021-01-01")
range(ModObs$time)
[1] "2012-01-01" "2020-12-01"
ModelErrorPlot_YF85 <- ggplot(data = ModObs, aes(x = time)) +
geom_ribbon(aes(ymin = low, ymax = high, fill = pond_name)) + # Adjust fill colors
geom_line(aes(y = med), color = "black") + # Black lines for the median values
facet_wrap(~ pond_name) + # Facet by pond_name
scale_fill_manual(values = color_mapping) + # Apply custom color palette for ribbons
labs(
x = "Year",
y = "Water Temperature (°C)",
color = "Pond Name",
fill = "Pond Name") + # Labels for legend
theme_classic() +
theme(plot.title = element_text(hjust = 0.5, size = 18),
axis.title = element_text(size = 22),
axis.text = element_text(size = 20),
axis.text.x = element_text(angle = 65, hjust = 1),
legend.position = "none",
strip.text = element_text(size = 16))
ModelErrorPlot_YF85
rsession-arm64(14805) MallocStackLogging: can't turn off malloc stack logging because it was not enabled.
ggsave("Corr_Error_YF85.png", plot = ModelErrorPlot_YF85, width = 9, height = 6)
Clear the environment
rm(list = ls())